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Abstract—Life-history studies of spe- 
cies discarded in fisheries are a low 
priority, particularly those of age and 
growth estimation. Consequently, almost 
everything is unknown about such spe- 
cies despite their having been caught 
as bycatch over the long term. In this 
study, age estimates were obtained by 
using the vertebrae of round stingrays 
(Urobatis halleri). We fitted 3 differ- 
ent growth models (von Bertalanffy, 
Gompertz, and logistic) to length-at-age 
data. Bayesian estimation of the various 
growth parameters was done by using 
the Markov chain Monte Carlo algo- 
rithm. Prior distributions of the param- 
eters disc width at birth (DW,) and the 
theoretical maximum disc width (DW...) 
were included and considered infor- 
mative. The priors for the growth coef- 
ficient (k) and the completion growth 
parameter (k.) were set as noninforma- 
tive (Rk for the von Bertalanffy growth 
function and k, for the Gompertz and 
logistic growth models). Results of our 
analyses indicate that growth bands 
are annual and that the round sting- 
ray lives up to 8 years. According to the 
Watanabe—Akaike information criterion 
for model selection, the von Bertalanffy 
growth function for sexes combined was 
selected as the best model. The mean 
values of marginal posteriors were as 
follows: DW,=9.52 cm (95% credible 
interval [CI]: 9.25—9.79), DW.,=32.50 cm 
(95% CI: 30.60-34.46), and k=0.114 
year | (95% CI: 0.101-0.129). 


Manuscript submitted 28 January 2022. 
Manuscript accepted 9 June 2022. 

Fish. Bull. 120:205—217 (2022). 

Online publication date: 8 July 2022. 
doi: 10.7755/FB.120.3-4.2 


The views and opinions expressed or 
implied in this article are those of the 


author (or authors) and do not necessarily 


reflect the position of the National 
Marine Fisheries Service, NOAA. 


Bayesian estimation of the age and growth 
of the round stingray (Urobatis halleri) in the 
Gulf of California in Mexico 


J. Fernando Marquez-Farias' 
L. Daniel Carrillo-Colin? 
Allan Rosales-Valencia' 


Raul E. Lara-Mendoza? 
Oscar G. Zamora-Garcia (contact 
author)* 


Email address for contact author: ozamora@sirbaa.com 


' Facultad de Ciencias del Mar 
Universidad Aut6noma de Sinaloa 
Paseo Claussen s/n 
Colonia Los Pinos 
82000 Mazatlan, Sinaloa, Mexico 


? Posgrado en Ciencias del Mar y Limnologia 
Universidad Nacional Autonoma de Mexico 
Avenida Ciudad Universitaria 3000 
04510 Coyoacan, Ciudad de Mexico, Mexico 


3 Direcci6n General Adjunta de Investigacién 
Pesquera en el Atlantico 
Instituto Nacional de Pesca y Acuacultura 
Secretaria de Agricultura y Desarrollo Rural 
Avenida Mexico 190 
Colonia Del Carmen 
04100 Coyoacan, Ciudad de Mexico, Mexico 


* Servicios Integrales de Recursos Bioldgicos, 
Acuaticos y Ambientales 
Avenida Macarela 3302 
Edificio Cardones |-104 
Palmilla Grand Residencial 
82112, Mazatlan, Sinaloa, Mexico 


The geography and oceanographic 
characteristics of the Gulf of California 
(GOC) provide a highly productive area 
that supports some of the largest fish- 
eries in Mexico. As a result, the GOC 
is one of the most productive regions 
for elasmobranchs. According to offi- 
cial fisheries statistics, in 2018, ~59% 
and ~39% of the landings of sharks and 
rays (batoids), respectively, came from 
the GOC (CONAPESCA’). 

Both the industrial and artisanal 
fisheries of the GOC contribute to the 
landings of batoids (Marquez-Farias, 
2002; Bizzarro et al., 2009). In addition 
to the current exploitation by industrial 


! CONAPESCA (Comisién Nacional de Acua- 


cultura y Pesca). 2018. Anuario estadistico 
de acuacultura y pesca 2018, 277 p. Comis- 
idn Nacional de Acuacultura y Pesca, 
Mazatlan, Mexico. [Available from website.] 


and artisanal ray fisheries, elasmo- 
branch populations have been affected 
by the shrimp-trawl fishery that has 
been operating in the GOC since 1921 
(Lluch-Cota et al., 2007). 

The most important species in the 
landings of the ray fishery belong to 
the families Rhinobatidae, Dasyati- 
dae, Myliobatidae, and Gymnuridae 
(Saldana-Ruiz et al., 2016). In addi- 
tion, because the fishing grounds of the 
artisanal ray fishery are inshore, the 
landings frequently include juveniles 
and pregnant females of many species 
(Marquez-Farias and Blanco-Parra, 
2006; Bizzarro et al., 2009). 

The composition of ray species in 
landings varies depending on the type 
of fishery, essentially because of the sea- 
sonal migrations of individuals, the hab- 
itat they occupy, and their susceptibility 
to capture. Given this situation and 
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the magnitude of the exploitation of elasmobranchs, the 
research priorities in the Mexican national plan of action 
for sharks focus on the more abundant species in the catch 
(CONAPESCA-INP, 2004). Considering the time series of 
historical catch of rays, the mixed-species fisheries, and the 
disadvantages of the low biological productivity of elasmo- 
branchs (Musick, 1999), the negative effect of multiple fish- 
eries on ray populations is presumably high. 

A recent study of species composition and fishing 
dynamics revealed the effect of the industrial shrimp 
fishery on the elasmobranch populations of the GOC 
(Garcés-Garcia et al., 2020). Results from that study 
confirm the importance of American round rays (family 
Urotrygonidae) as bycatch. The round stingray (Urobatis 
halleri) is the most abundant fish species in the bycatch 
of the industrial shrimp fishery in the GOC (L6pez- 
Martinez et al., 2010; Rabago-Quiroz et al., 2011). Unfor- 
tunately, up to 95% of the bycatch volume of the shrimp 
fishery in the GOC is discarded at sea because of its low 
commercial importance (Zamora-Garcia, 2015), limiting 
the study of the life history of the round stingray and 
other species. 

The round stingray is distributed from California (in 
the United States) to Ecuador. It is a benthic and small 
(maximum disc width [DW]: 30.8 cm) batoid inhabiting 
soft mud or sand bottoms and is usually found in near- 
shore waters less than 15 m deep but can occur at depths 
down to 91 m (McEachran and Notarbartolo-di-Sciara, 
1995). Aggregations are modulated by water temperature 
(>17°C) (Ebert, 2003; Hoisington and Lowe, 2005). Its 
food preferences include benthic invertebrates, bivalves, 
and small fishes. It is a viviparous species, and brood size 
ranges from 1 to 6 pups per season, depending on the size 
of the female (Babel, 1967), after 3 months of gestation 
(Nordell, 1994). 

Differences in life-history traits in some elasmobranch 
species, such as the shovelnose guitarfish (Pseudobatos 
productus) and Pacific angel shark (Squatina californica), 
have been documented between populations in California 
and those farther south on the western coast of Baja Cal- 
ifornia and in the GOC in Mexico (Marquez-Farias, 2007, 
2020). The growth characteristics of round stingrays have 
been estimated for only the California population and 
are unknown for the GOC population (Hale and Lowe, 
2008). The round stingray is categorized as a species of 
least concern in the IUCN Red List of Threatened Species 
(Lyons et al., 2015). From a regional perspective, popula- 
tions of the round stingray and other rays likely have been 
affected by the shrimp-trawl fishery, given its intensity for 
more than 70 years in the GOC. 

In this study, we report the results of Bayesian analy- 
ses used to investigate the age and growth of the round 
stingray in the GOC for the first time. First, using stan- 
dard aging protocols, we used sectioned vertebrae to 
interpret and count the growth bands and to test the 
accuracy of band counts by readers. Then, we fitted 3 non- 
linear growth models to age—length data pairs to esti- 
mate growth. Selection of the model that best explained 
and predicted the length-at-age data was based on the 


Watanabe—Akaike information criterion (WAIC), an 
appropriate approach for Markov chain Monte Carlo 
numerical integration. Growth parameters are presented 
as posterior probabilities that summarize the informa- 
tion of data and the prior probability of the parameters. 
The results of this study provide insight into key popula- 
tion characteristics necessary to model the demography 
of the round stingray and to test its potential resilience 
to fishing mortality. 


Material and methods 
Data and sample collection 


Specimens of the round stingray were collected monthly 
during 2007-2012 as bycatch from trawl surveys con- 
ducted by the artisanal shrimp fishery in Teacapan, 
Sinaloa, Mexico (22°32’N, 105°45’W). We recorded the sex 
of each specimen and measured the DW (in centimeters) 
in a straight line, avoiding the curvature of the body. In 
addition, a section of 8—10 cervical vertebrae was extracted 
from each specimen to prevent intravertebral variability 
of band pairs (Natanson et al., 2018). Finally, the samples 
were labeled and frozen at —20°C at lab facilities for sub- 
sequent analysis. 


Sample processing 


Excess tissue and neural arches were removed from 
the vertebral segments manually. Each centrum was 
removed and mounted onto a flat piece of wood with 
synthetic resin. The vertebral diameters (VDs) were 
measured sagittally (to the nearest 0.05 mm) by using 
a Wixey” WR100 digital vernier (Barry Wixey Develop- 
ment, Seattle, WA). Using an IsoMet Low-Speed Saw 
(Buehler Ltd., Lake Bluff, IL), with twin blades sepa- 
rated by 0.24 mm, we cut each vertebra in a sagittal 
plane through the focus. Each vertebral section was 
mounted on a slide and photographed by using trans- 
mitted light with an Olympus SZ61 stereoscopic micro- 
scope (Olympus Corp., Tokyo, Japan) equipped with an 
Infinity1-2 digital camera (Teledyne Lumenera, Ottawa, 
Canada). We created an image bank for the interpreta- 
tion and counting of growth bands. 


Size frequency and the relationship of size to vertebral 
diameter 


We prepared DW-frequency histograms for comparison 
of size between the sexes. We used a t-test to determine 
whether there was a difference in the mean size between 
the male and female specimens. The relationship between 
VD and DW was assessed by using linear regression. 


2 Mention of trade names or commercial companies is for identi- 
fication purposes only and does not imply endorsement by the 
National Marine Fisheries Service, NOAA. 
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Age estimation 


In the vertebral sections of round stingrays, growth- 
band pairs appear opaque and translucent under trans- 
mitted light (Cailliet et al., 2006). Age estimates were 
made by counting band pairs along the corpus calcareum 
after identifying the birthmark, which can be seen as a 
change in the angle close to the focus. After establishing 
the aging criteria, 3 independent readers counted the 
growth bands. 


Aging bias and precision 


Reproducibility of counts of growth bands was assessed 
by using the index of average percent error (IAPE) (Equa- 
tion 1) and the coefficient of variation (CV) (Equation 2). 
Next, the accuracy was evaluated by using the precision 
index (D) (Equation 3) for estimates from the 3 readers 
together and between readers (Chang, 1982; Goldman 
and Musick, 2006). Finally, the consensus of growth-band 
counts was visually evaluated by using age-bias plots 
(Campana et al., 1995). The following equations were 
used: 


IAPE = 100| -. zy 


ast(istbcsll 


sand (2) 


o-feel sh : 


where N = the number of individuals aged; 
R = the number of readers; 
x;; = the age estimate z for individual j; 
x; = the mean age calculated for individual j; 
A = the number of agreements; and 
B = the number of readings done. 


For further analyses, we considered only age estimates 
on which at least 2 readers agreed; otherwise, the read- 
ing was discarded (Goldman et al., 2012). After growth 
bands were counted, the radius (distance from the focus 
to the edge) of a vertebra and each band’s thickness were 
recorded by using Image-Pro Plus image analysis software 
(vers. 6.0, Media Cybernetics Inc., Rockville, MD). 


Periodicity of band formation 


The periodicity of growth-band formation in round sting- 
rays was assessed through centrum edge analysis (CEA) 
and marginal increment analysis (MIA) following Cailliet 
(2015). First, we categorized the edge type of each section 
of a vertebra as opaque or translucent. Then, the relative 


frequency of each edge type was tabulated in bimesters 
(2-month intervals). The MIA (Equation 4) was con- 
ducted bimonthly by using the ages 2-8 years, following 
Natanson et al. (1995). Bimesters with mean marginal 
increment (MI) values close to 1 were interpreted to be 
the time of year when the growth cycle is about to be com- 
pleted (Cailliet et al., 2006): | 


MI = At a } (4) 


Tn ~Tp-1 


where VR = the vertebral radius; 
r, = the distance from the focus to the last band 
pair; and 
r,-1 = the distance from the focus to the penultimate 
band pair. 


One-way analysis of variance (ANOVA) was used to assess 
both the strength of the difference in total band counts 
between each pair of readers and the differences in the MI 
among bimesters (Gerrodette, 2011; Carrillo-Colin et al., 
2021). 

A Bayesian approach was used for the t-test, the linear 
regression, and the ANOVA described above. The Bayes 
factor (BF',,.) was used to measure the weight of the evi- 
dence of the alternative hypothesis (H,) against the null 
hypothesis (H,) (Kass and Raftery, 1995; Rouder et al., 
2009): 


BE ge Hil Eg: (5) 


The null hypothesis corresponds to the effect size (5) of 0, 
and the alternative hypothesis was specified as a Cauchy 
prior with location and scale parameters of 0.000 and 
0.707, respectively (Jeffreys, 1961). 

All statistical treatments were done under the Bayesian 
approach by using the programs JAGS (vers. 4.3.0; 
Plummer, 2003) and JASP (vers. 0.15; JASP Team, 2021), 
statistical software R, vers. 4.0.5 (R Core Team, 2021), and 
the BayesFactor R package (vers. 0.9.12-4.2; Morey and 
Rouder, 2018). We relied on the guidelines proposed by 
Jeffreys (1961) for interpreting BF). 


Growth modeling 


To maximize the accuracy of the length-at-age predic- 
tions for both sexes of round stingrays and to contrast 
different hypotheses about growth, we used the von Ber- 
talanffy, Gompertz, and logistic growth models (Katsa- 
nevakis and Maravelias, 2008). Growth models included 
in this study were those versions that incorporate disc 
width at birth (DW,), as recommended for rays (Smart 
et al., 2016): 

von Bertalanffy growth function (VBGF), 


DW, = DW., - (DW., - DW, Je“; (6) 
Gompertz growth model, 


DW, = DW, ;and (7) 
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logistic growth model, 


kot 
ie 


Be, (8) 
"DW. + DWyek* 


where DW, = the estimated disc width (in centimeters) at 
a given age ¢ (in years); 
DW.,,=the theoretical maximum disc width (in 
centimeters); 
k = the growth coefficient (year”*); and 
k, = the completion growth parameter (year™'). 


Model fitting was done by using the Bayesian approach, 
which requires an initial probability distribution (prior) 
for each model parameter. We conducted a meta-analysis 
based on growth parameters published in the literature 
for round stingrays (Table 1). The mean and standard 
deviation of published growth parameters were used to 
set the parameters of the prior distributions (informa- 
tive). The priors for the parameters k and k, were chosen 
to be relatively noninformative as the uniform distribu- 
tion, with a minimum of 0 and maximum of 10 to avoid 
negative values and let the model acquire the needed 
curvature (Smart and Grammer, 2021). The prior for 
DW, was informative and set as a normal distribution 
(ND): the population mean equal to 33.05 (standard 
deviation [SD] 10.07). The prior for DW, was also infor- 
mative with an ND (population mean=7.32, o=0.29). The 
prior of the residual o was uniform, with a minimum of 
0 and a maximum of 100. 

The posterior distribution for the parameters of the 
3 growth models (Equations 6-8) was computed by 
using the Bayes theorem (Gelman et al., 2003). With 
the actual growth parameters substituted into the 
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Bayes formula and the terms simplified, the equation 
is as follows: 


P(DW.,,k,,, DWp|d) 


co 9 "nN ? 


_ P(d|DW,,,k,,DW,)P(DW.,,k, ,DWo) 


co 9 "nN ? 


~ [ P(d|DW.,, ky, DWo)P(DW., kz, ;DWo) 38 


co 9 'N ? Ce) e102) 


(9) 


where P = the probability; 
k,, = the growth parameter of the n model; 
d = the data; and 
00 = the partial derivate for each parameter 0. 


The complexity of this multiparameter equation requires 
numerical integration. Therefore, we used a Markov 
chain Monte Carlo algorithm to obtain the posterior 
probability distribution by drawing random samples 
from the vector of all parameters of each model. For this 
purpose, we used JAGS, assuming an additive error. To 
assess the posterior distribution convergence, 4 Markov 
chains were simulated, with 100,000 iterations, a burn-in 
of 5000 runs, and a thinning of 1 iteration (Carrillo-Colin 
et al., 2021). Trace plots were assessed visually (Gelman 
and Rubin, 1992). 


Model selection 


To decide which model best explained and predicted the 
growth of round stingrays with Markov chain Monte Carlo 
numerical integration, we used the WAIC. This Bayesian 
criterion considers the entire parameter’s likelihood 
matrix (Watanabe, 2009). First, model selection was made 
after comparing the WAIC values of the 3 competing mod- 
els, and the one with the lowest WAIC value was selected. 
After that, the selected model, including the version with 


Table 1 


Reported growth parameters for male (M) and female (F) round stingray (Urobatis halleri) or for both sexes com- 
bined from other studies. These published means for disc width at birth (DW,), theoretical maximum disc width 
(DW.,), and the growth coefficient (Rk) were used in this study to set the parameters of the prior distributions and 
to compare the von Bertalanffy growth function (VBGF) and the Gompertz growth model for round stingray in the 
southern Gulf of California in Mexico. A single asterisk (*) indicates a value of DW,, determined through the trans- 
formation of the theoretical maximum total length (TL,,.) of 47.2 cm with the equation DW=1.33+0.53(TL) obtained 
in this study. A double asterisk (**) indicates values for the completion growth parameter from the Gompertz 
growth model. A dash indicates that the k value was not provided by the author. 


DW.. DW, k 


Model (cm) (cm) (year') 
VBGF 31.00 7.50 
25.00 7.50 
VBGF Peis (a 47? 
VBGF 22 AD 7.60 
28.56 7.00 
Gompertz 46.56 7.00 
30.75 7.00 


Location of study 
Seal Beach, California 


Gulf of California, Mexico 
Seal Beach, California 


Southern Gulf of 
California, Mexico 


Source 
Babel (1967) 


Morales-Azpeitia et al. (2013) 
Hale and Lowe (2008) 


Diliegros-Valencia (2019) 
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sex variation, was fitted, and the WAIC was com- 


puted. The WAIC values of the selected model for. 


both sexes (reduced model) and its version that 
includes sex variation (full model) were then com- 
pared to decide if there was evidence of differ- 
ences in growth between sexes. 


Results 


In total, 244 round stingrays were examined 
(99 females and 145 males). We found weak 
evidence (BF,,)=1.00) in the mean differences 
in DW between sexes, where males were 1.01 
cm (SD 0.50) larger than females (95% credible 
interval [CI]: 0.03—2.01). The size range for sexes 
combined was 7.8—24.5 cm DW (mean: 14.36 cm 
DW [SD 3.87]) (Fig. 1). We found conclusive evi- 
dence (BF ,)>100) of a positive linear correlation 
between DW and VD (coefficient of determina- 


€ 
= 
eo 
oS 
= 
2) 
2 
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2.0 ES 3.0 
Vertebral diameter (mm) 


Figure 2 


Relationship between vertebral diameter and disc width for round 
stingrays (Urobatis halleri) sampled during 2007—2012 in the south- 
ern Gulf of California in Mexico (for sexes combined). Dashed lines 
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tion [r7]=0.87) with regression coefficients a and 
b of 18.43 (95% CI: 12.04—24.74) and 55.26 (95% 
CI: 52.62—58.03), respectively, described by the 
equation DW=18.43+55.26(VD) for sexes com- 
bined (Fig. 2). 

Growth-band counts ranged from 0 to 8. The 3 readers 
agreed on 38% of the counts. When just 2 readers agreed 
(62% of cases), their count differed by 1 band in 59% of 
cases and by 2 bands in 3% of cases. The indicators of 
bias and precision for the total sample were as follows: 
IAPE=20%, CV=27%, and D=56%. Visual inspection of 
the age-bias plots between each pair of readers revealed 
minimal variation around the 1:1 ratio (Fig. 3, A-C). The 
results of the post hoc tests done as part of the ANOVA 
provide anecdotal evidence that the mean counts between 


Frequency 
o a 


Disc width (cm) 


Figure 1 


Length—frequency distributions for female (white bars) and male 
(black bars) round stingrays (Urobatis halleri) sampled during 
2007-2012 in the southern Gulf of California in Mexico (number of 
rays=244). 


indicate the 95% credible interval. 


readers 1 and 2 were different (BF',)=2.02). In contrast, for 
the mean counts between readers 1 and 3 (BF’;,=0.27) and 
readers 2 and 3 (BF,,=0.17), the evidence of difference was 
null (Fig. 3D). 

Results from the CEA indicate that the highest fre- 
quency of opaque edges occurred in the summer months, 
likely the time when deposition of growth bands began. 
Results from the MIA indicate that growth-band depo- 
sition appears close to being completed in November— 
December when the maximum MI was reached (Fig. 4). 
The lower MI estimated in the fourth bimester 
(July-August) was probably caused by the low 
sample size. However, we found substantial evi- 
dence in the mean differences in MI between 
bimester 2 (March—April) and bimester 4 (July— 
August), bimester 3 (May—June) and bimester 
4 (July-August), and bimester 5 (September— 
October) and bimester 6 (November—December) 
(Suppl. Table) (online only). Results from both 
CEA and MIA indicate that every growth band 
was formed annually. 

The estimated ages for sexes combined ranged 
from 0 to 8 years. The mean age of all the sam- 
ples was 2.27 years (SD 1.87), with 23.8% and 
17.6% of individuals aged at 0 and 1 years, 
respectively. As it is with many elasmobranchs, 
overlap in size at age increased with age for our 
specimens of the round stingray, which can be 
seen in the size-at-age key in Table 2. 

The model selected according to the WAIC as 
the best fit to the data was the VBGF (Table 3). 
When sex was included in the VBGF, the WAIC 
increased by 6.46 units (from 804.47 to 810.93). 
Given the loss of parsimony by including sex as 
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Figure 3 


Age-bias plots of counts of growth bands on vertebrae of round stingrays (Urobatis halleri) sampled between 
2007 and 2012 in the southern Gulf of California in Mexico. Black circles represent mean counts (A) of 
reader 2 (R2) relative to reader 1 (R1), (B) of reader 3 (R3) relative to reader 1, (C) of reader 3 relative to 
reader 2, and (D) between readers 1-3 from ad hoc Bayesian one-way analysis of variance. Error bars repre- 
sent the range of band counts assigned by readers. The dashed line in panels A-C indicate 1:1 equivalence. 
Reproducibility and precision of age estimates were assessed by using the index of average percent error 
(IAPE), coefficient of variation (CV), and precision index (D). 


a factor, we selected the VBGF fitted with data for the 
sexes combined. The mean values of marginal posteriors 
from the VBGF were as follows: DW)=9.52 cm (95% CI: 
9.25-9.79), DW.,=32.50 cm (95% CI: 30.60—-34.46), and 
k=0.114 year’ (95% CI: 0.101-0.129). The growth curve 
is monotonic as expected. However, because the sample 
size decreased with age, the model lost accuracy at ages 6 
and 8 (Fig. 5). The VBGF for the round stingray predicts 
growth from the birth size of 9.52 cm DW to 11.92 cm DW 
for age 1 and to 14.21 cm DW for age 2. Therefore, the 


growth rate from age 0 to age 1 and age 2 was 2.48 cm/year 
and 2.21 cm/year, respectively. 

The marginal posterior probability distribution of the 
VBGF parameters (Fig. 6) indicates a smooth and well- 
defined shape for every parameter (Fig. 6, A, C, E). The 
relationship of DW.,, to k is typically inverse (Fig. 6B). Sim- 
ilarly, the relationship of k to DW, is inverse but with a 
slope that is less steep (Fig. 6F), whereas the relationship 
of DW, to DW.,, is dispersed and centered without correla- 
tion (Fig. 6D). 
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Figure 4 


The periodicity of band formation was assessed by using marginal 
increment (MI) analysis and edge percentage analysis for round 
stingrays (Urobatis halleri) sampled during 2007-2012 in the 
southern Gulf of California in Mexico. The black circles represent 
the posterior mean marginal increments, and the error bars indi- 
cate the 95% credible intervals. The gray columns represent the 
percentages of opaque edges. Numerals above the circles indicate 
number of rays. Results of analyses are shown for each bimester 
(2-month interval): January—February, March—April, May—June, 
July-August, September—October, and November—December. 


Table 2 


Discussion 


The maximum age of the round stingray found 
in this study (8 years) is the same as that found 
by Babel (1967) but lower than that reported by 
Hale and Lowe (2008) (14 years). Furthermore, 
the latter authors found statistical differences in 
growth between males and females. This differ- 
ence could be related to the size structure of both 
samples. Although the size range obtained in this 
study is similar to that in the study by Hale and 
Lowe (2008), their sample was composed mainly 
of males (71%), and our sample did not have such 
bias in sex proportion (59% males). Differences 
in availability could produce the phenomenon 
of apparent change in the growth rate like gear 
selectivity can (Walker et al., 1998). Also, the dif- 
ference in conclusions between the studies could 
be a result of the use of the different statistical 
approach in our study (Gerrodette, 2011). 
Compared to other members of the family 
Urotrygonidae, the round stingray has a max- 
imum age that is similar to that of the Pana- 
mic stingray (Urotrygon aspidura) from the 
Pacific coast of Columbia (Torres-Palacios et al., 
2019) but is shorter than that reported for the 


Age-size key for round stingrays (Urobatis halleri) sampled between 2007 and 2012 in the southern Gulf of Cal- 
ifornia in Mexico (for sexes combined). In the center of the table are the numbers of individuals that correspond 
to each disc width (DW) class and age class. In the margin of the table, the total number of stingrays (n) and the 
mean and standard deviation (SD) for each DW class and age class are presented. 


Age (years) 


A 


9-10 
10-11 
11-12 
12-13 
13-14 
14-15 
15-16 
16-17 
17-18 
18-19 
19-20 
20-21 
21-22 
22-23 
24—25 


n 58 43 37 36 35 23 
Mean (cm) 9.7 ER RY 14.0 16.6 18.5 19.4 
SD 0.89 1.19 1.55 1.32 0.95 1.20 


11 
19.4 
LAG 


Mean 
(years) 


0.0 
0.0 
0.1 
0.3 
1.0 
1.3 
1.6 
2.0 
2.4 
2.8 
3.4 
4.3 
4.7 
4.9 
5.2 
5.8 
8.0 
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Table 3 


Comparison of the von Bertalanffy, Gompertz, and logistic growth models fit to length-at-age 
data for round stingrays (Urobatis halleri) sampled between 2007 and 2012 in the southern Gulf 
of California in Mexico. The parameters are disc width at birth (DW,, in centimeters), theoret- 
ical maximum length (DW.,., in centimeters), and the growth coefficient (R, year“) for the von 
Bertalanffy growth function and the completion growth coefficient (Ro, year‘) for the Gompertz 
and logistic growth models. Posterior values and 95% credible intervals (lower quartile: 2.5%; 
upper quartile: 97.5%) are provided from the fit of the growth models with a Markov chain Monte 
Carlo algorithm to data pooled for both sexes. The Watanabe—Akaike information criterion (WAIC) 
was used in model selection. The deviance information criterion (DIC) is reported to facilitate 
comparison of parameter estimates with those from previous studies in which the frequentist 
approach was used. SD=standard deviation; SK=standard error. 


Credible interval 


Model Parameter 2.5% 97.5% WAIC (SE) 


von Bertalanffy 804.47 (20.3) 
820.77 (21.1) 


Logistic 897.75 (21.5) 


blotched stingray (Urotrygon chilensis) 
(14 years) in the Mexican Pacific Ocean 
(Guzman Castellanos, 2015). This new 
information about the maximum age 
of the round stingray has considerable 
implications for studies of population 
dynamics. 

The evidence of a relationship between 
VD and DW supports the use of vertebrae 
for aging (Goldman et al., 2012). The pre- 
cision values (IAPE, CV, D) for samples in 
our study are within the expected inter- 
val for this group of species (Campana, 
2001) and agree with the readings that 
have been done for the family Urotrygo- 


Disc width (cm) 


nidae (Mejia-Falla et al., 2014; Guzman 4 5 
Castellanos, 2015; Santander Neto, 2015; noe. Waats) 
Diliegros-Valencia, 2019). Figure 5 

Results from both the CEA and MIA Estimated von Bertalanffy growth curve for round stingrays (Urobatis halleri) 
support the assumption of an annual in the southern Gulf of California in Mexico (for sexes combined). The black 
pattern of growth-band formation, and (male) and open (female) diamonds represent observed lengths at age for 
annual formation has also been reported stingrays sampled during 2007-2012. The gray shaded area indicates the 95% 


for round stingrays in California (Hale credible interval. 
and Lowe, 2008) and in the southeastern 
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Figure 6 


The Bayes posterior probability distributions of the (A) theoretical maximum disc 
width (DW.,), (C) disc width at birth (DW,), and (E) growth coefficient (2), the param- 
eters of the von Bertalanffy growth function, paired with scattergrams of the rela- 
tionships between the model’s parameters, (B) DW., to k, (D) DW, to DW.,, and (F) k 
to DW,. Values are predicted from the 100,000 Markov chain Monte Carlo simula- 
tions used to fit the von Bertalanffy growth function to observed length-at-age data for 
round stingrays (Urobatis halleri) sampled during 2007-2012 in the southern Gulf of 


California in Mexico. 


GOC (Diliegros-Valencia, 2019). The opaque bands were 
deposited in the summer months. They are broad and rep- 
resent rapid growth, and the translucent bands formed 
in winter are thinner and reflect slower growth (Cailliet 
et al., 1986; Cailliet and Goldman, 2004). Evidence of vari- 
ability between individuals in the number of band pairs 
has been documented for the round stingray (James and 
Natanson, 2020), and such a source of bias was prevented 
in this study by using just the cervical vertebrae for all 
round stingrays. Also, band-pair counts without validation 
need to be taken with caution. Although there are more 
direct methods for validating the periodicity of growth 
bands of aquatic organisms, CEA and MIA are advanta- 
geous because they yield reliable and comparable results 
(Campana, 2001; Carrillo-Colin et al., 2021). Although 
new methods and approaches will rise to overcome the 
current controversies in age studies for batoids (James 


and Natanson, 2020), our statistical approach and results 
set a baseline for the age and growth of the round stingray 
in the tropical region of the Pacific Ocean. 

In this study, we explored 3 conceptually different mod- 
els that include DW,, an appropriate parameter for vivip- 
arous chondrichthyans (Cailliet et al., 2006). We avoided 
exploring other models with more parameters because 
the WAIC penalizes them or because they do not have 
any biological significance (Richards, 1959; Katsaneva- 
kis, 2006). The continuous growth pattern observed in 
this study (as indicated by parameter estimates from the 
VBGF) and the representativeness of most age groups 
indicate that the round stingray completes its life cycle in 
the study area (Bergstad et al., 2021). Interestingly, the 
round stingray has been reported to segregate by size and 
life stage (Silva-Garay and Lowe, 2021), but there are no 
signs of segregation in our data set. Size segregation and 


214 


length selectivity of fishing gears can bias growth data 
to indicate a discontinuous pattern (Walker et al., 1998) 
that can favor more flexible models (with an inflection 
point) or even 2-phase growth models arguing investment 
in reproductive processes (Torres-Palacios et al., 2019). 

The round stingray is not commercially targeted in the 
GOC artisanal fisheries (Bizzarro et al., 2009). However, 
the low sample size for large individuals (>6 years and 
>22 cm DW) might indicate an effect of fishing mortality 
on its population due to being caught as bycatch of shrimp 
trawlers (Garcés-Garcia et al., 2020), as implied by the 
difference between the observed maximum DW (24.50 cm) 
and the estimated mean DW,, (32.50 cm). Fortunately, the 
statistical framework used in this study allowed us to esti- 
mate DW. by combining the prior information and sup- 
port of data considering the parameter’s uncertainty in a 
probabilistic manner. 

Size selectivity does not seem to be an explanation for 
the absence of large individuals because the maximum 
DW obtained during this study is very similar to those 
of other studies (27.00 cm; Hale and Lowe, 2008), and 
such a minimal difference in maximum DW does not 
indicate a significant difference in swimming speeds or 
turning ability. The latter may be true for species with 
higher swimming ability and size, such as the golden 
cownose ray (Rhinoptera steindachneri), the eagle rays 
(Myliobatis spp.), or those of the genus Hypanus, which 
are more likely to outswim the relatively slow shrimp 
trawls. In contrast, the poor mobility of round stingrays 
limits their ability to escape, and their flattened shape 
restricts size selection through the trawl. Consequently, 
fishing mortality is likely to affect the size range of the 
round stingrays available in the fishing area and sus- 
ceptible to being captured by the trawlers. Although our 
surveys were conducted with small trawl nets used in 
the artisanal shrimp fishery, the size range sampled was 
convenient for this study and included individuals in 
all classes for both sexes. Furthermore, similar sizes for 
the round stingray were recently reported from fishery- 
independent surveys conducted with large vessels and in 
different areas, depths, and seasons in the GOC (Garcés- 
Garcia et al., 2020). 

In this study, we investigated the round stingray’s age 
and growth and developed a fully Bayesian approach 
for estimation of all quantities. In contrast, leaving the 
model fitting to the data alone could produce misleading 
estimates of parameters, especially when the observed 
data do not follow the expected pattern, are incomplete, 
or are biased. The stochastic approach in which Monte 
Carlo simulation is used can be advantageous for rem- 
edying the lack of fit. However, it also requires assign- 
ing a probability distribution to the parameters without 
weighting them with the data to produce the final dis- 
tributions of the parameters (Neer et al., 2005). Alter- 
natively, the Bayesian approach in this study relies on 
the explicit use of uncertainty for the growth parame- 
ters. A parameter’s uncertainty is expressed as a prob- 
ability distribution and combined with the information 
contained in the data (likelihood) throughout the Bayes 
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theorem (Equation 9); it is possible to obtain an enriched 
marginal posterior probability distribution (Fig. 6) of each 
parameter (Gelman et al., 2003). This fact is essential in 
age and growth studies when parameters can be informa- 
tive based on the knowledge of the life history of a species, 
as is the case for DW, (Carrillo-Colin et al., 2021; Smart 
and Grammer, 2021). 

Neither Monte Carlo simulation nor the Bayesian 
approach can correct for biases in the data caused by 
length-selective fishing mortality. Troynikov and Walker 
(1999) stated that the influence of length-selective mor- 
tality on determining growth parameters is not easily 
explained because it is a combination of 2 overlapping 
effects on lengths for different age classes. These authors 
proved that a growth pattern can be affected by an appar- 
ent change in growth rate due to length-selective fishing 
mortality. 

Reporting the parameters as probability distributions 
(and their credible intervals) is not only fundamentally 
different from classical (frequentist) estimation but also 
useful as inputs in further quantitative demographic 
analysis or stock assessment for elasmobranchs (Cortés, — 
2002). However, the philosophy of measuring the statis- 
tical support of each growth parameter when fitting the 
model goes beyond the scope of our work described here 
(Gerrodette, 2011). 

In our study, we used the Bayes factor in the t-test and 
ANOVA to evaluate the size of the effects, and we used 
the WAIC for the selection of the best growth models. The 
WAIC is convenient for more complex stochastic numeri- 
cal approximations. 

The convergence of the growth parameters DW..,, k, and 
DW, was satisfactorily achieved on the basis of their mar- 
ginal posterior distributions (Carrillo-Colin et al., 2021). 
The informative prior distribution for DW, favored model 
fitting because it helped in the selection of realistic birth 
size values (Rolim et al., 2020; Smart and Grammer, 2021). 
The k value (from the VBGF) we obtained in our study 
ranged between 0.11 and 0.13 year™’ and is close to the 
range of 0.2-0.5 year’ that has been reported for species 
of the families Rhinobatidae, Torpedinidae, and Urotrygo- 
nidae (Cailliet and Goldman, 2004). However, it must be 
noted that this broad difference in k values can also be 
related to sample size, aging method, and verification, val- 
idation, and growth model fitting techniques (Cailliet and 
Goldman, 2004). 


Conclusions 


In the GOC, chondrichthyan populations have long been 
affected by fisheries that trawl for teleost and shrimp 
species. Capture as bycatch in those fisheries causes 
significant mortality for species of no commercial value, 
including for the round stingray and other demersal 
small rays (Garcés-Garcia et al., 2020). Variations in 
life-history characteristics among populations can be 
expected in chondrichthyan species because of the fol- 
lowing: 1) cold-water species tend to grow larger and 
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probably live longer than tropical water species; 2) dif- 
ferences in food availability, habitat type, and environ- 
mental conditions affect their growth rate; and sample 
distribution can be biased 3) by the effects of migration 
and segregation and 4) by the length selectivity of fishing 
gears (Walker, 1998; Walker et al., 1998). In this study, 
we estimated growth parameters by using the counts of 
growth bands in vertebrae taken from sampled round 
stingrays. We found a maximum age of 6 or 8 years, a 
remarkably short life span for a chondrichthyan species, 
indicating that the round stingray has a short intrinsic 
generation time and, therefore, a high intrinsic rate of 
population increase (King and McFarlane, 2003). Stock 
assessment and ecological risk analyses require inputs 
from age, growth, and reproductive rate studies. The 
Bayesian estimation of this study can be used in future 
quantitative analyses to gain the advantage of consid- 
ering uncertainty in the estimation of parameters (Punt 
and Hilborn, 1997). 


Resumen 


Los estudios sobre la historia de vida de las especies 
descartadas en pesquerias son de baja prioridad, espe- 
cialmente aquellos sobre estimaciones de edad y crec- 
imiento. En consecuencia, se desconoce casi todo sobre 
estas especies, a pesar de haber sido capturadas de 
forma incidental por mucho tiempo. En este estudio, las 
estimaciones de edad se obtuvieron utilizando las vérte- 
bras de la raya redonda (Urobatis halleri). Se ajustaron 
3 diferentes modelos de crecimiento (von Bertalanffy, 
Gompertz y logistico) a datos de longitud a la edad. La 
estimacién Bayesiana de los distintos parametros de 
crecimiento se realiz6 mediante el algoritmo cadenas 
de Markov de Monte Carlo. Se incluyeron y consideraron 
distribuciones previas informativas de los parametros 
ancho del disco al nacer (DW,) y ancho del disco maximo 
tedrico (DW.,,). Las distribuciones previas para el coefi- 
ciente de crecimiento (2) y el parametro de terminacién 
del crecimiento (&,) se establecieron como no informa- 
tivas (k para la funcién de crecimiento de von Berta- 
lanffy y k, para los modelos de crecimiento Gompertz 
y logistico). Los resultados de nuestros an4lisis indican 
que las bandas de crecimiento son anuales y que la raya 
redonda vive hasta 8 afios. Segtn el criterio de infor- 
macioén de Watanabe—Akaike para la seleccién del mod- 
elo, la funcién de crecimiento de von Bertalanffy para 
sexos combinados fue seleccionada como el mejor mod- 
elo. Los valores promedio de las posteriores marginales 
fueron los siguientes: DW,)=9.52 cm (intervalo de credi- 
bilidad [IC] al 95%: 9.25—-9.79), DW,.=32.50 cm (IC 95%: 
30.60-34.46), y k =0.114 afio' (IC 95%: 0.101—-0.129). 
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Abstract—Management of recreational 
fishing for greater amberjack (Seriola 
dumerilt) in the Gulf of Mexico involves 
size regulations and closed seasons. 
Water temperature, salinity, fight and 
handling times, and barotrauma can 
influence survival of released fish. We 
examined postrelease mortality and 
behavior by using acoustic telemetry 
of the movements of 78 fish, moni- 
tored for up to 58 d across 3 sampling 
efforts in 2018, 2019, and 2020 over a 
depth gradient of 29-64 m. To assess 
descender devices as a mitigation tool, 
we assigned fish to 2 treatments: sur- 
face release without swim bladder 
venting and release with a descender 
device. Cox proportional hazards mod- 
els were used to assess the effects of 
site depth, release treatment, bait type 
(jigging or live bait), fishing injury, tag- 
ging injury, fight and handling times, 
surface and bottom temperatures, and 
fish length. We found neither a positive 
association between mortality risk and 
site depth, as might be expected from 
barotrauma, nor increased survivor- 
ship for fish released with a descender 
device. The best-supported model con- 
sidered only fish length as a factor in 
postrelease mortality; legal-size fish 
(>864 mm in fork length) had a mortal- 
ity risk 20 times greater than that of 
smaller fish. Our results indicate that 
sublegal-size fish released because of 
size restrictions face much lower mor- 
tality risk than legal-size fish. 
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The greater amberjack (Seriola dumerilt) 
is a cosmopolitan tropical and temper- 
ate predatory fish taken in recreational 
and commercial fisheries (Smith-Vaniz, 
2002). In the Gulf of Mexico, the greater 
amberjack is considered overfished 
and currently experiencing overfish- 
ing (SEDAR, 2020). Past amendments 
to the Gulf of Mexico Fishery Man- 
agement Council fishery management 
plan included requirements for the use 
of non-stainless-steel circle hooks with 
natural baits, an increase in the rec- 
reational minimum size from 762 mm 
(30 in) to 864 mm (84 in) in fork length 
(FL), and seasonal closures during open 
recreational seasons for red snapper 
(Lutjanus campechanus) and other reef 
fish species (SEDAR, 2020). These fac- 
tors make it likely that more greater 
amberjack will be caught at sublegal 
sizes (<864 mm FL at the time of the 
study) and out of season and then 
released. Therefore, there is a need to 
gain a better understanding of factors 
that contribute to discard mortality. 
Greater amberjack discarded in rec- 
reational fisheries face several potential 


sources of mortality. These sources 
include the following: 1) at-vessel mor- 
tality (AVM), meaning a death after 
a fish is hooked but before landing; 
2) capture and handling mortality 
(CHM), meaning a lethal event during 
capture and handling after landing that 
prevents successful release; and 3) post- 
release mortality (PRM). 

Several factors have been predicted 
to affect discard mortality associated 
with recreational fishing for greater 
amberjack. Swim bladder barotrauma 
is associated with increased PRM risk 
for some reef fishes (Curtis et al., 2015; 
Runde and Buckel, 2018). Swim blad- 
ders in acanthomorph fishes, like the 
greater amberjack, lack a pneumatic 
connection to the gut and are termed 
physoclistous (Lagler et al., 1962). 
Physoclistous fishes cannot expel gas 
from the swim bladder lumen through 
the pharynx and as a result must resorb 
gas into the blood stream through 
the oval (Woodland, 1913). Therefore, 
depth-related pressure changes from 
rapid fishing ascents may cause injury 
to swim bladder tissues, increase stress, 
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and make descent to deeper habitats difficult for fish after 
release. Results from a previous study of red snapper indi- 
cate that barotrauma is associated with both immediate 
discard mortality (AVM in our study) and delayed mortal- 
ity (PRM in our study) (up to 72 h postrelease) and that 
barotrauma can be mitigated by venting the swim bladder 
or using a descender device to recompress the swim bladder 
and then releasing the fish at depth (Curtis et al., 2015). 

Other factors may contribute to AVM and PRM in reef 
fishes like the greater amberjack. Fight time during fish- 
ing may raise risk of AVM and PRM by increasing stress 
(Mohan et al., 2020), depleting energy, and elevating preda- 
tion risk. In addition, fight time may increase with fish size 
and fishing depth, and the combined factors may contribute 
to mortality. Water temperature, salinity, and dissolved oxy- 
gen may also be additional sources of physiological stress. 

Biotelemetry was used to assess PRM in greater amber- 
jack at artificial reefs off the coast of Alabama (Jackson 
et al., 2018). Using depth-logging acoustic tags to study 
fish that were vented and released, Jackson et al. (2018) 
found a PRM estimate (18.8%) similar to the estimate for 
1 of 3 modeled PRM scenarios in a stock assessment from 
the same time (SEDAR, 2014). In our study, we caught 
greater amberjack in the northern Gulf of Mexico using 
common recreational angling techniques in 2018, 2019, 
and 2020. We used direct observation to calculate AVM 
and CHM and used data from acceleration- and depth- 
logging acoustic tags to estimate PRM of fish released at 
the surface or with a descender device. The main study 
objectives were 1) to estimate survivorship and associated 
PRM of discarded greater amberjack and 2) to evaluate the 
following factors as potential predictors of PRM risk: fish 
size, fight and handling times, use of a descender device, 
site depth, capture method (live bait or jigging), observed 
fishing and tagging injuries, and water temperature at the 
surface and on the bottom. 


Materials and methods 
Receiver deployment 


Vemco VR2AR! acoustic receivers (Innovasea Systems 
Inc., Boston, MA) were deployed at 16 sites on artificial 
reefs before each study period in 2018, 2019, and 2020. 
Artificial reefs were steel and concrete pyramids, sunken 
boats, a barge, an oil rig jacket, a fuel tank, a grain hopper, 
and a submarine (Suppl. Table 1) (online only). One or 2 
receivers were placed near each site (Suppl. Table 1) (online 
only). For deployment, VR2AR receivers were attached to 
cement moorings with approximately 2 m of polypropy!- 
ene line attached to the detachable lug, and a midline 
swivel was present between the mooring and lug to pre- 
vent twisting. Two non-compressible 20-cm (8-in) trawl 
floats were attached to the receiver collar with 3.5-m of 


1 Mention of trade names or commercial companies is for identi- 
fication purposes only and does not imply endorsement by the 
National Marine Fisheries Service, NOAA. 


polypropylene line with a midline swivel to provide buoy- 
ancy. Each assembled receiver, along with its moorings and 
floats, was deployed at the sea surface and released with 
the manufacturer’s acoustic release function with a Vemco 
VHTx-69k transponding hydrophone (Innovasea Systems 
Inc.) and VR100 receiver (Innovasea Systems Inc.) at the 
end of each study period. 


Fish collection and tagging 


Tagging took place on 16-17 and 23-24 August 2018, on 
30 April and 13-14 May 2019, and on 17-18 and 21-22 
August 2020. These dates were chosen on the basis of the 
availability of personnel and vessels and the aim to pro- 
vide an opportunity to examine mortality under different 
water temperatures and at different fish sizes, which var- 
ied among the 3 study periods. 

Greater amberjack were caught and tagged by following 
common recreational fishing methods on chartered boats 
(FV Lady Ann and FV Escape). Fish were caught by using 
2 methods opportunistically, either with live bait (11/0 or 
12/0 circle hooks, 340—450-¢g weights, and a 36.3—45.4-kg 
test monofilament leader with a swivel tied to the main- 
line) or with artificial lures (140-—170-g jig heads with soft 
plastic lures). A data recorder on board, using a stopwatch, 
noted the time from when an angler first had a fish on 
the line to when the fish was landed (fight time, in sec- 
onds) and the time from the landing to the release of a 
fish (handling time, in seconds). Upon their landings, fish 
were immediately measured and ventilated with a salt- 
water hose. Standard, fork, natural total, and stretch total 
lengths were recorded to the nearest millimeter. Injuries 
to fish were noted prior to release, and all dead fish prior 
to landing (AVM events) were recorded. 

Fishing release treatment (surface or descender device) 
was alternated in the order that fish were landed. Before fish- 
ing began each day, we randomly chose the starting release 
treatment. For the release treatment in which a descender 
device was used, a SeaQualizer Descending Device (stan- 
dard model, SeaQualizer, Davie, FL) was attached to the 
fish’s lower jaw and set to the deepest possible release 
depth (increments of 15.2 m, 30.5 m, and 45.7 m) at each 
site, requiring a release depth less than the site depth. The 
fish and descender device were lowered with approximately 
2.3 kg of weight attached to a fishing rod along with a dig- 
ital video camera (Hero5 Black Edition, GoPro Inc., San 
Mateo, CA) for observing potential predation. 

Results from a previous study (Jackson et al., 2018) 
indicate that the observed “release condition” of a fish 
was indicative of PRM risk. We assigned a number for the 
release condition (Patterson et al., 2001; Jackson et al., 
2018) of fish released at the surface as follows: 1, fish 
immediately oriented and swam downward rapidly; 2, 
fish appeared disorientated and swam down slowly; 3, fish 
appeared disoriented and remained at the surface for sev- 
eral minutes; and 4, fish was dead and unresponsive at the 
surface. It was not possible to assess release condition for 
fish released with a descender device because fish drifted 
out of the view of the camera in most cases. 
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Fish were tagged with Vemco V9AP depth and accel- 
erometry acoustic transmitters (Innovasea Systems Inc.) 
coded with the following parameters: activity algorithm, 
triaxial root mean square of acceleration, 5 samples/s, 
20—40-s delay, accelerometer range of +4.9 m/s, high power 
setting, a slope (resolution) of 0.3 m for depth, depth range 
of 68 m, and acceleration-to-depth transmission ratio of 
1:1. We initially attached transmitters externally on fish 
rather than intraperitoneally to avoid inadvertent swim 
bladder venting that could result from surgical implanta- 
tion and that might interfere with testing of barotrauma 
effects. In 2018, acoustic transmitters were attached exter- 
nally on a second external dart tag (Floy FH-69 stainless- 
steel dart tag, Floy Tag Inc., Seattle, WA), with the tag 
attached to the dart with a 68-kg monofilament line. Dart 
tags were inserted and locked to interneurals (pterygio- 
phores) below the spiny dorsal fin on the left side. Prior to 
tagging, acoustic transmitters were attached to tags with 
marine epoxy and 2 zip ties, with the transmitter’s longi- 
tudinal axis parallel to the external tag. A second dart tag 
(FIM-96 nylon dart tag, Floy Tag Inc.), with a unique tag 
ID, phone number, and website URL on it, was placed cau- 
dally to the stainless-steel dart tag on the fish’s left side. 

Because of acoustic transmitter shedding (see the 
“Results” section) in 2018, we modified the tagging pro- 
cedure for 2019 and 2020. Vemco V9AP accelerometer 
transmitters were placed intraperitoneally through a 
small incision, just wide enough to pass the transmitter. 
The incision was made in the left abdomen, just dorsal 
(1-2 cm) and anterior (1—2 cm) to the cloaca. The swim 
bladder was not visible in any tagging surgeries, and no 
evidence of venting was observed. After tag implanta- 
tion, incisions were closed with 2 interrupted sutures by 
using monofilament suture thread. In addition to acoustic 
transmitters, fish received a dart tag with ID and contact 
information as described above: FIM-96 tags were used in 
2019, and FH-69 tags were used in 2020. 


Postrelease fate inferred from acoustic telemetry 


We used depth and acceleration data from tagging of fish 
to infer their fate (Whitney et al., 2016). We tabulated raw 
depth and acceleration data from receivers by time and 
examined graphs of all raw data. We plotted raw data over 
an approximately 5-d period (if available) at the end of 
the detection period for each fish or, in the cases of sus- 
pected mortality or tags that had been shed (see the “Post- 
release fate based on acoustic telemetry” subsection in the 
“Results” section), over the time period of suspected death 
or tag shedding. Fish were determined to be dead when 
transmissions indicated that the tag was on the bottom 
and acceleration had ceased. Rapid onset of high acceler- 
ation values were predicted to indicate possible predation 
or carcass scavenging. 


Survivorship 


We estimated overall survivorship, accounting for all sources 
of discard mortality (AVM, CHM, and PRM); at-vessel 


survivorship, accounting for AVM and CHM; and postre- 
lease survivorship, accounting for PRM. Survivorship (S) 
and its standard error (SE) were estimated from acoustic 
telemetry (see the “Results” section) and data on recapture of 
fish by using equations 17 and 18 of Pollock and Pine (2007): 


= (1) 
n 


SE(S) = y= (2) 
Nv 


where x = the number of surviving fish; and 
n = the number of fish: in the context of this equa- 
tion, total number of fish caught or dead prior 
to landing for estimates of overall survivorship, 
number of fish landed for at-vessel survivor- 
ship, or number of fish released for postrelease 
survivorship. 


We produced Cox proportional hazards models (Cox, 
1972) to test for the contribution of factors (i.e., fishing 
related variables, abiotic conditions, release treatment, and 
fish size) to postrelease mortality. Cox models were pro- 
duced by using the survival package (vers. 3.2-13; Thernau, 
2021) in R (vers. 4.1.1; R Core Team, 2021). For all models, 
we report the statistic from the likelihood ratio test in which 
models were compared to a null model and the P-values for 
each test and their associated model predictors. Separate 
models were run for externally tagged fish and internally 
tagged fish because many fish in the former group shed 
their transmitters. Only some fish were recaptured, and 
because it is not possible to determine death events after 
acoustic tags ceased to transmit data for fish that were not 
recaptured, we restricted the observation period of our Cox 
models and survivorship curves to the period when acous- 
tic tags were transmitting and receivers were deployed. 
This restricted observation period did not change model 
outcomes but limited the duration of survivorship curves 
to <60 d. All 9 recaptured fish had transmitters that were 
detected until the end of predicted transmitter battery life 
or receiver retrieval, but 4 of these fish were externally 
tagged fish that shed their acoustic tags. Acceleration and 
depth data from the tags attached to these 4 fish indicate 
that tags were shed within range of receivers while they 
were still transmitting (see the “Results” section). 

We initially considered the following potential predic- 
tors: FL (in millimeters); fight and handling times; release 
treatment (categorical, fish released at surface or with 
descender device); site depth; bait type (categorical, live 
bait or jigging); fishing injury (categorical, no injury versus 
visible injury from fishing gear before handling and tag- 
ging); tagging injury (categorical, no injury versus visible 
bleeding from site of external tagging); surface, mid-depth, 
and bottom water temperatures; surface, mid-depth, and 
bottom salinities; and surface, mid-depth, and bottom dis- 
solved oxygen. Barotrauma injuries (e.g., expanded swim 
bladders and bulging eyes) were not included in models 
because they were never observed, although internal tis- 
sue damage may have been present. Release condition 


Boyle et al.: Evaluation of postrelease mortality of Seriola dumenii with depth and acceleration data from tags 221 


was not tested because that information was not available 
for fish released with a descender device. Before includ- 
ing predictors in Cox models, we tested for correlation 
among predictors with Spearman rank correlation tests 
(Kneebone et al., 2021) (Suppl. Table 2) (online only). Many 
abiotic variables (e.g., temperature, salinity, and dissolved 
oxygen) were highly correlated (Suppl. Table 2) (online only); 
therefore, we reduced the analysis to include bottom and 
surface temperature as the only abiotic variables. Several 
additional correlations were observed among other predic- 
tors included in initial Cox models (e.g., fight time and fish 
length; Suppl. Table 2 [online only]), and these correlations 
are considered in our interpretation of results. 

The internal placement of tags in fish in 2019 and 2020 
may have caused additional handling stress, and although 
attempts were made to not vent the swim bladder, we 
inadvertently may have partially vented it. Therefore, we 
used separate Cox models to analyze fish with transmit- 
ters placed externally (in 2018) and those with transmit- 
ters placed intraperitoneally (in 2019 and 2020). Stepwise 
model selection was conducted with the stepAIC function 
from the package MASS (vers. 7.3-54; Venables and Ripley, 
2002) in R, and both forward and backward selections were 
attempted to select the model with the lowest Akaike infor- 
mation criterion (AIC). We used the AIC with correction for 
small sample sizes (AICc) (Burnham and Anderson, 2002) 
and considered models with a difference in AICc from that 
of the best model (AAICc) of 2 units or less to have compara- 
ble support, unless the model with more variables differed 
by addition of a single parameter, indicating a potentially 
spurious variable (Burnham and Anderson, 2002). 

Predictor data were missing for several fish: fight time 
(1 fish from 2018), handling time (1 fish from 2018), bottom 
temperature (14 fish from 2019 and 2020), and fight and 
handling times (6 fish from 2019 and 2020). Therefore, to 
explore all factors, we began with full models for the subset 
of fish with available predictor data and completed step- 
wise AICc model selection. In one instance, for internally 
tagged fish with bottom temperature data available, the 
full model was overparameterized and did not converge. We 
ran exploratory models with single predictors and removed 
predictors with the highest AlCc before running the step- 
wise procedure on the fullest model that would converge. 
We then ran separate procedures for stepwise AIlCc model 
selection on the complete data sets (externally tagged fish 
from 2018 [number of fish (n)=23] and internally tagged 
fish from 2019 and 2020 [n=55]) that did not include the 
predictors that were missing for some fish. 

Kaplan—Meier survival function curves were generated 
for variables informed by Cox models by using the sur- 
vfit function in the survival package in R and were plotted 
with the survminer package (vers. 0.4.9; Kassambra et al., 
2021) in R. 


Results 


We tagged 78 greater amberjack at 15 artificial reef sites 
that ranged in depth from 29 to 64 m during separate 


efforts in 2018, 2019, and 2020 (Table 1). Three fish died 
before successful release (Table 1). Fish D1 died on the 
boat and was bleeding severely from the hook site (classi- 
fied as a CHM). Fish D2 died from a propeller injury, and 
fish D3 died from predation by a shark before landing 
(both considered an AVM). All other fish released at the 
surface had a release condition of 1. In underwater video 
from the descender rig, fish were usually out of view 
because of a long leader between the camera and fish. 
Therefore, fish released with a descender device near the 
bottom were rarely observed, and it was not possible to 
determine their release condition. Further, no predation 
or potential predators (e.g., carcharhinid sharks) were 
observed on video. 

Fourteen acoustically tagged fish were detected on 
receivers at reefs where they were not originally tagged 
(Table 1). Three of these 14 fish, fish 118, 120, and 121, 
were detected on 2 reefs at close proximity (0.5 km apart) 
throughout the study. Six fish tagged in 2020 (fish 58, 58, 
104, 105, 107, 115) were not detected on the first day after 
being tagged and released (Table 1). 


Postrelease fate based on acoustic telemetry 


Patterns of acceleration and depth over time were dis- 
tinct among fish determined to have died versus lived. 
Seven fish acoustically tagged during 2018—2020 were 
inferred to have died (Fig. 1, Table 2). Profiles from tags 
of dead fish indicate only brief depth changes and accel- 
eration oscillations greater than 2 m/s” (Fig. 1, Table 2). 
Fish 22, 31, 51, 117, and 119 were dead very soon after 
release (<2 h), and fish 23 and 47 remained alive longer 
(Fig. 1, Table 2). In contrast, fish inferred to have lived 
had strong and consistent depth and acceleration oscil- 
lations until the end of predicted transmitter battery life 
or deployments of acoustic receiver arrays, until they 
emigrated from receiver arrays, or in the case of some 
fish with externally placed transmitters (as was done in 
2018), until their tag was shed (Fig. 2, Table 2, Suppl. 
Figs. 1 and 2, Suppl. Table 2 [online only]). 

Acceleration and depth profiles indicate that 13 of 23 
fish tagged in 2018 shed their tags, and information on 
the recapture of fish by anglers confirms this outcome in 
4 instances (Table 1); secondary external tags were not 
shed. Results of analysis of the carcass of a recaptured fish 
indicate that the stainless-steel dart remained in place 
but that the monofilament line attachment point broke on 
the dart, perhaps as a result of drag on the tag or the fish 
actively scraping the tag off. No tags that were intraperi- 
toneally placed in 2019 or 2020 were shed. Telemetry pro- 
files for fish inferred to have shed tags indicate relatively 
consistent depth and acceleration oscillations that abruptly 
ended at the time of the presumed tag shedding (Fig. 2, 
Table 2, Suppl. Figs. 1 and 2, Suppl. Table 3 [online only]). 
The tag records for fish inferred to have shed tags differed 
from the tag records of fish presumed to have died; for fish 
that died, acceleration and depth variation was evident but 
less consistent during the brief period before their appar- 
ent death (Fig. 1). Ten fish were detected briefly (<4 h) 
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Table 1 


Summary of tagging effort and the fork lengths (F'Ls) and fates for greater amberjack (Seriola dumerili) caught on charter boats in 
the northern Gulf of Mexico during 2018—2020. Tagged fish were released either at the surface (S) or with a descender device (DC). 
In 2018, fish had an acoustic tag attached to a primary external tag (Ex.), and in 2019 and 2020, fish had an acoustic transmitter 
placed intraperitoneally (In.), in addition to a secondary external tag. Symbols after the last detection date indicate if the tag was 
shed (st) or if the fish was recaptured (7). Fish <864 mm FL were of sublegal size. Either jigging (J) or live bait (L) was used to 
catch fish. Fates of fish include alive (A), dead (D), lived until recaptured (R), emigrated outside of detection range during the life of 
the transmitter battery and before receivers were retrieved (E), moved between sites with receivers (M), and unknown because the 
fish spent <68 h in the receiver array (U). Symbols in the “Fate” column indicate that the external tag was shed, either determined 
from acceleration and depth data (*) or confirmed from angler recapture (**). Fish D1 was dead before release, fish D2 died from a 
propeller injury before landing, and fish D3 died of predation before landing. UK=unknown because no data were taken. 


Site Detection date 
Fish Release Tag depth Tagging ee er Ce FP FL Bait FT HT 
ID type type (m) date First Last (mm) type (s) (s) Fate 
01 S Ex. 30.5 16-Aug-2018 16-Aug-2018 23-Aug-2018" 950 L 99 146 Ny 
02 DC Ex. 33.2 16-Aug-2018 16-Aug-2018 16-Aug-2018 1006 i 172). 126 U 
03 S Ex. 33.2 16-Aug-2018 16-Aug-2018 1-Aug-2019' 750 He FB). wash Oi R- 
04 DC Ex. 37.0 16-Aug-2018 16-Aug-2018 5-Jun-2019" 995 Li 250 155 R- 
05 S Ex. 37.0 16-Aug-2018 16-Aug-2018  17-Sep-2018* 630 J 145 135 Ne 
06 DC Ex. 37.0 16-Aug-2018 16-Aug-2018  12-Sep-2018* 527 J 54 UK A 
07 S Ex. 37.0 16-Aug-2018 16-Aug-2018 23-Aug-2018 1100 L UK 29 A 
08 DC Ex. 45.0 16-Aug-2018 16-Aug-2018 21-Aug-2018" 895 D 141 159 A 
09 S Ex. 45.0 16-Aug-2018 16-Aug-2018 16-Aug-2018 930 1 108 72 i 
10 DC Ex. 45.0 16-Aug-2018 16-Aug-2018 19-Sep-2018" 940 L 173 107 Rie 
11 S Ex. 45.0 16-Aug-2018 16-Aug-2018 16-Aug-2018" 900 i 171 58 A 
12 DC Ex. 49.0 16-Aug-2018 16-Aug-2018 8-Oct-2018 855 IL 117 Ae A 
13 S Ex. 52.4 17-Aug-2018 17-Aug-2018 17-Aug-2018 984 I: 134 98 AE 
14 DC Ex. 52.4 17-Aug-2018 17-Aug-2018 17-Aug-2018 1029 ib 185 75 U 
15 S Ex. 52.4 17-Aug-2018 17-Aug-2018 17-Aug-2018 986 i 167 93 U 
16 DC Ex. 38.0 17-Aug-2018 17-Aug-2018 17-Aug-2018 896 iB £70! | 07 AE 
17 S Ex. 38.0 17-Aug-2018 17-Aug-2018 17-Aug-2018 917 ji 269 62 U 
18 DC Ex. 38.0 17-Aug-2018 17-Aug-2018 17-Aug-2018 905 L 99 210 U 
19 S Ex. 38.0 17-Aug-2018 17-Aug-2018 17-Aug-2018 966 J 260 58 AE 
20 DC Ex. 38.0 17-Aug-2018 17-Aug-2018 17-Aug-2018" 919 L 154 93 A 
1 S Ex. 33.2 23-Aug-2018 23-Aug-2018 5-Jun-2019" 836 L 98 129 Be 
9 DC Ex. 33.2 23-Aug-2018 23-Aug-2018 23-Aug-2018 779 ii 100 157 D 
23 S Ex. 37.0 23-Aug-2018 23-Aug-2018 26-Aug-2018 540 J 55 70 D 
24 DC In. 29.1 30-Apr-2019 30-Apr-2019 25-Jun-2019 835 J UK UK A 
25 S In. 29.1 30-Apr-2019 30-Apr-2019 23-May-2019 736 L UK UK AE 
26 DC In. 29.1 30-Apr-2019 30-Apr-2019 6-Jan-2020' 730 Ts UK R 
27 S In. 29.1 30-Apr-2019 30-Apr-2019 15-Oct-2021' 730 Ie UK R 
28 DC In. 64.0 30-Apr-2019 30-Apr-2019 2-Aug-2019' 829 iy UK R 
Zo S In. 64.0 30-Apr-2019 30-Apr-2019 23-May-2019 825 L UK AE 
30 DC In. 64.0 30-Apr-2019 1030 L 156 318 U 
31 S In. 64.0 30-Apr-2019 30-Apr-2019 30-Apr-2019 983 im 374 278 D 
32 DC In. 32.3 13-May-2019 13-May-2019 25-Jun-2019 718 J 99 320 AM 
33 S In. 32.3 13-May-2019 13-May-2019 25-Jun-2019 735 L 10 219 A 
34 DC In. 32.3 13-May-2019 13-May-2019 31-May-2019 707 | 112 - "998 AE 
35 S In. 32.3. 13-May-2019 13-May-2019 25-Jun-2019 720 L 90 136 AM 
36 DC In. 36.9 13-May-2019 13-May-2019 25-Jun-2019 634 J 90 285 A 
37 S In. 36.9 13-May-2019 13-May-2019 25-Jun-2019 703 L 60 189 A 
38 DC In. 64.0 14-May-2019 14-May-2019 22-Jun-2019 650 J 96 534 AE 
39 S In. 64.0 14-May-2019 14-May-2019 22-Jun-2019 809 if 115 355 AE 
40 DC In. 64.0 14-May-2019 14-May-2019 24-May-2020! 705 9 120 430 R 
Al S In. 64.0 14-May-2019 14-May-2019 25-Jun-2019 800 J 165 394 A 
42 DC In. 64.0 14-May-2019 14-May-2019 22-Jun-2019 734 i 108 538 AE 
43 S In. 29.1 14-May-2019 14-May-2019 25-Jun-2019 656 L AT 248 A 
4A DC In. 29.1 14-May-2019 14-May-2019 27-May-2019 716 J 100 445 AE 
A5 S In. 29.1 14-May-2019 14-May-2019 27-May-2019 705 J 130 548 A 
AT DC In. 37.0 17-Aug-2020 17-Aug-2020 27-Aug-2020 961 is 177. 924 D 


(Continued on next page) 
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Release Tag 


In. 
In. 
In. 
In. 
In. 
In. 
in: 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
In. 
Tr 
In. 
In. 
In. 
In. 


Site 


depth 


(m) 


37.0 
37.0 
36.1 
36.1 
34.8 
34.8 
34.8 
34.8 
34.8 
34.8 
34.8 
34.8 
34.8 
34.8 
34.8 
34.8 
34.8 
38.0 
38.0 
38.0 
38.0 
38.0 
38.0 
38.0 
38.0 
36.1 
36.1 
36.1 
47.8 
47.8 
47.8 
47.8 
37.0 
37.0 
33.2 


Tagging 
date 


17-Aug-2020 


17-Aug-2020 
17-Aug-2020 
17-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
21-Aug-2020 
21-Aug-2020 
21-Aug-2020 
22-Aug-2020 
22-Aug-2020 
22-Aug-2020 
22-Aug-2020 
16-Aug-2018 
23-Aug-2018 
16-Aug-2018 


Table 1 (continued) 


Detection date 


First 


17-Aug-2020 
17-Aug-2020 
17-Aug-2020 
17-Aug-2020 
18-Aug-2020 
25-Aug-2020 


18-Aug-2020 
30-Aug-2020 
18-Aug-2020 
28-Aug-2020 
18-Aug-2020 


28-Aug-2020 
30-Aug-2020 
18-Aug-2020 
29-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 
18-Aug-2020 


18-Aug-2020 
18-Aug-2020 
22-Aug-2020 
21-Aug-2020 
21-Aug-2020 
22-Aug-2020 
22-Aug-2020 
22-Aug-2020 
22-Aug-2020 


Last 


18-Aug-2020 
14-Sep-2020 
21-Aug-2020 
17-Aug-2020 
1-Oct-2020 
23-Sep-2020 


18-Aug-2020 
1-Oct-2020 
1-Oct-2020 
1-Oct-2020 
1-Oct-2020 


14-Sep-2020 
1-Oct-2020 
1-Oct-2020 
26-Sep-2020 
1-Oct-2020 
1-Oct-2020 
18-Aug-2020 
15-Sep-2020 


26-Aug-2020 
1-Oct-2020 
24-Aug-2020 
14-Sep-2020 
21-Aug-2020 
25-Sep-2020 
22-Aug-2020 
1-Oct-2020 
22-Aug-2020 
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(Fig. 3). Few detections occurred for these fish, but the 
depth and acceleration values were similar to fish that 
lived (Fig. 3, Table 2). Five fish were never detected by any 
receivers (Table 1). 

Among released fish that died, 4 greater amberjack 
(57%) were released at the surface and 3 fish (48%) were 
released with a descender device. For externally tagged 
fish for which PRM was not observed, the length of time 
during which their tags were detected in 2018 was vari- 
able (median: 7.8 d; 1st quartile: 2.8 d; 3rd quartile: 33.0 d) 
because some fish shed their transmitters and were not 
recaptured and some likely emigrated away from reefs 
with acoustic receivers. For internally tagged fish, the 
length of time during which their tags were detected in 
2019 and 2020 were longer (median: 38.9 d; 1st quartile: 
14.1 d; 3rd quartile: 43.8 d). Eight fish were captured 90 d 
or later (maximum: 899 d) after tagging. 


Recapture rate and postrelease survivorship 


The fish recapture rate from recreational and commer- 
cial angling (fishing mortality rate) was 11.5%. In 2018, 
acoustic tags were placed externally. Postrelease mortality 
was inferred from acceleration and depth data for 2 fish in 
2018, 1 fish in 2019, and 4 fish in 2020, or 9% of fish tagged 
in all 3 years combined. 

We estimated survivorship with data for all fish in the 
study except fish that either were never recaptured or 
were detected acoustically for less than 68 h, the latest 
observed PRM event (Fig. 1). Overall survivorship after 
accounting for all sources of discard mortality, AVM, 
CHM, and PRM (i.e., including the 3 fish that died before 
release), was moderate (S=84.6% [SE 4.1]). At-vessel sur- 
vivorship following AVM and CHM was 94.8% (SE 4.4), 
and postrelease survivorship was 88.7% (SE 4.2). 
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Figure 1 


Acceleration and depth profiles from acoustic telemetry for tagged greater amberjack (Seriola dumerili) in 
the northern Gulf of Mexico inferred to have died after release. Table 2 provides the rationale for the fate 
and time of inferred death of each fish. Fish were caught and tagged on charter boats during 2018-2020. 


Cox proportional hazards models: 
factors of postrelease mortality 


The 2 best models used to examine PRM of fish tagged 
with external transmitters in 2018 received comparable 
support (AAICc: 0.86), and both included bait type and 
surface temperature as predictors (Table 3). In the second- 
best model, bait type and surface temperature were signif- 
icant predictors (Table 3); however, the very small hazard 
ratios associated with these predictors indicate a tiny 
decrease in mortality associated with fish caught on a jig 
and an even smaller decrease in mortality with each 1°C 


increase in surface temperature. No predictors were sig- 
nificant in the best model, which also included fishing 
injury as a variable (Table 3). Effects of fight and handling 
times were tested by using a subset of data from 2018 for 
fish for which these values were available, and these fac- 
tors were not significant predictors of mortality in any of 
the models used in this study (Suppl. Tables 4—6) (online 
only). The Kaplan—Meier survivorship curve for fish tagged 
externally in 2018 (Fig. 4) indicates that survival probabil- 
ity dropped to 94.7% 2 h postrelease and to 88.8% 68 h 
postrelease and then remained at that level throughout 
the rest of the period of acoustic monitoring. 
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Table 2 


Fates of 25 greater amberjack (Seriola dumerili) tagged and released between 2018 and 2020 in the northern Gulf of Mexico, 
inferred with depth and acceleration data from acoustic tags. Fates include alive (A), alive and detected on 2 separate reef sites 
before emigrating away from the receiver array (AME), alive and emigrated (AE), alive and shed tag (AST), alive but detected 
for <68 h (U), alive and detected on 2 separate reef sites for 68 h (UM), and dead (D). Also noted for each fish are the figure (Fig.) 
showing acceleration and depth data that were recorded before and continued to be recorded until the last acoustic detection or 
fate assignment and either the time from release to fate assignment determined from acoustic telemetry data (alive, dead, or shed 
tag) or the date the fish was recaptured. 


Time to 
fate (h) or 
Fish recapture 
ID Fate Rationale for fate assignment Fig. date 
22 D Brief acceleration activity and little depth variation indicates fish remained near bottom 1 2.00 
(at depth of ~36.5 m) and ventured up to a depth of only 30 m. Appears to have died 
within 2 h postrelease. Movement seen in slightly deeper water 20 h after first detection 
consistent with movement of carcass by scavengers or of tag by currents. 
23 D Appears near bottom from 34 to 67 h postrelease, but acceleration indicates fish still moving 1 67.97 
vigorously. Makes one last foray to shallow water (depths of 38.5—15.5 m) ~67 h postre- 
lease (perhaps dragged by a predator and dropped). Tag remains on bottom after 68 h. 

31 D Brief, strong acceleration quickly ceases. Depth values drop quickly, indicating that fish is 1 0.28 

dead on the bottom. 

47 D Acceleration stops suddenly at time of death with additional brief peaks (~20, ~45, and 1 7.62 

~70 h postrelease) not associated with depth changes but consistent with movement of 
carcass by scavengers or of tag by currents. 

51 D Acceleration stops suddenly, and fish is at the bottom (depth of ~39 m) within 1.20 h 1 1.20 

postrelease. 
117 D Detected on the bottom (depth of ~40 m), and acceleration stopped within 0.20 h postre- d 0.20 
lease. Additional, isolated accelerations at ~20, 30, and 94 h postrelease consistent with 
movement of carcass by scavengers or of tag by currents. 
119 D Appears dead within 1 min of release. Only one acceleration transmission indicative of fish 1 0.02 
movement and first depth detections occur on bottom at depth of ~51.5 m. Additional 
movement ~6 h postrelease is consistent with movement of carcass by scavengers or of 
tag by currents. 
O01 AST ‘Tag shed (inferred) 166 h postrelease. Consistent variation of acceleration and depth and 2 166.38 
an immediate cessation of movement when tag is on bottom, at a depth of ~33 m. Note 
similarity to depth and acceleration profiles of fish 03, 08, and 10 in Figure 2. 
03 AST Tag shed, confirmed from recapture by angler. Tag appears to have been shed 164 h 2 1-Aug-2019 
postrelease. Note similarity to fates of fish 01, 08, and 10. 
08 AST Tag shed (inferred) 116 h postrelease. Note similarity to fates of fish 01, 03, and 10. 2 115.95 
10 AST  Tagshed, confirmed from recapture by angler during acoustic monitoring period. Tag appears 2 19-Sept-2018 
to have been shed 86 h postrelease. Note similarity to fates of fish 01, 03, and 08 in Figure 2. 
Data was right-censored at 800.50 h postrelease in Cox proportional hazards models. 

BF A Remains alive and emigrates, confirmed by angler recapture. Note the large gaps in detection 2 15-Oct-2021 

of tag transmissions at ~765—770 h, ~815—820 h, and ~840—850 h postrelease. 

29 AE Appears alive until it leaves acoustic receiver array. P 543.17 
118 A Appears alive through the end of tag transmissions. 2 813.00 
120 A Appears alive through the end of tag transmissions. 2 960.98 

02 U Very few detections, but fish appears alive when last detected. 3 0.60 

09 ie! Only 3 acceleration and 2 depth detections, but fish appears to have been alive when last 3 3.07 

detected. 

14 U Very few detections, but fish appears to have been alive when last detected. 3 1.72 

15 U Very few detections, but fish appears to have been alive when last detected. 3 1.22 

18 U Only 3 acceleration and 5 depth detections, but fish appears to have been alive when last 3 0.28 

detected. 

48 UM _ Few detections. Fish is detected at a second site after emigrating from first site. Fish 3 11.70 

appears alive when last detected. 

50 AME Fish is detected at a second site after emigrating from first site. Fish appears alive after 3 93.83 

leaving receiver array. 

55 U Very few detections, but fish appears to have been alive when last detected. 3 3.50 
115 U Very few detections, but fish appears to have been alive when last detected. 3 67.23 
121 U Very few detections, but fish appears to have been alive when last detected. 3 0.45 
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Figure 2 


Acceleration and depth profiles from acoustic telemetry for tagged greater amberjack (Seriola dumerili) 
in the northern Gulf of Mexico inferred to have lived after release until the end of the battery life for their 
transmitters, until they emigrated from the detection range of the acoustic receiver at their release site, or 
until they shed their tag. Table 2 provides the rationale for the fate of each fish. Fish were caught and tagged 


on charter boats during 2018-2020. 


For fish internally tagged in 2019 and 2020, fish length 
was the best predictor of survivorship, with PRM risk 
increasing with fish length (Table 4). The second-best 
model (AAICc: 0.67) included fishing injury in addition to 
fish length (Table 4). With legal length (2864 mm FL) 


considered a categorical predictor, the estimated hazard 
ratio indicates 20 times greater PRM risk for fish of legal 
lengths (Table 4, Fig. 5). Legal-size fish had an estimated 
survival probability of 55.6% 7.6 h postrelease, and 
sublegal-size fish had a survival probability of 97.6% from 
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Figure 3 


Acceleration and depth profiles from acoustic telemetry for tagged greater amberjack (Seriola dumerili) 
that were briefly detected with receiver arrays deployed at 16 artificial reef sites in the northern Gulf 
of Mexico. Fish were caught and tagged on charter boats during 2018-2020. All acceleration and depth 
data through the last detection of each fish are shown. The fish were inferred to be alive at the time of 
the last detection. All fish for which data are depicted in this figure, except fish 50, had a fate classified 
as unknown because their tags were detected for less than 68 h (the longest time observed in a case of 
postrelease mortality). 
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Table 3 


Results from analysis with the Cox proportional hazards models fit to data for a subset of greater amberjack (Seriola dumerili) 
tagged with external acoustic transmitters in 2018 in the northern Gulf of Mexico (number of fish [n]=23). In the full model, the 
following predictors (P) are considered: fork length (FL); use of a descender device or not (DC); bait type (BJ), live bait or jigging; 
fishing injury (FI), an injury attributed to fishing gear; tagging injury (TI), an injury associated with the tagging procedure; 
surface temperature (ST); site depth (D); and bottom temperature (BT). Stepwise model selection based on Akaike information 
criterion (AIC) scores was used to determine the model with the most support. For each model, the AIC corrected for small sam- 
ple sizes (AICc), the difference in AICc between the model and the best model (AAICc), and the likelihood ratio test statistic used 
to determine support of each model relative to a null model are provided. For each predictor in each model, the beta coefficient (B) 
and its standard error (SE) and the hazard ratio and its 95% confidence interval (95% CI) are provided. An asterisk (*) indicates 
that the predictor or model is significant (P<0.05). 


Likelihood 
Hazard Hazard ratio 
AlCe AAICe P ratio ratio 95% CI statistic df 


Full model 
26.29 19.02 
43.32 0.837 1x 10° 6x10" 


30,010 11,000 ~0 to oo 
3232 0.066 ~0 to oo 
17,330 2 x10" ~0 to o 
HOCIG. “oe xo" ~0 to oo 
1.6 x 10° 21.2 ~0 to 0 
23,560 0.7 x 10° ~0 to 0 
81,380 25 “107 ~0 to oo 


Second-best model 
8.13 0.86 
6.8x107° 4x10%1x10°% <0.001* 
5107" 32x 1081x102? F_.20,001" 


Best model 
7.26 0.00 
12,170 Oe“ 
52,850 Vet al {Ole 
25,950 hee 108 


Null model 
11.30 AA 


Survival probability > 
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Figure 4 


The Kaplan—Meier survivorship curve for greater amberjack (Seriola dumerili) with external 
acoustic tags in the northern Gulf of Mexico in 2018, (A) for the entire period (>50 d) for which 
transmitters and receivers were active and (B) for the first 100 h postrelease. Each cross along the 
curve indicates when data were censored because a fish emigrated away from acoustic receivers or 
shed its tag. The gray area indicates the 95% confidence interval for survival probability. 
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Table 4 


Results from analysis with the Cox proportional hazards models fit to data for greater amberjack (Seriola dumerili) 
implanted with internal transmitters in the northern Gulf of Mexico in 2019 and 2020 (number of fish [n]=55). In 
the models, the following predictors (P) are considered: fork length (FL); use of a descender device or not (DC); bait 
type (BJ), live bait or jigging; fishing injury (FI), an injury attributed to fishing gear; tagging injury (TI), an injury 
associated with the tagging procedure; site depth (D); surface temperature (ST); and fish length (LG), with fish 
categorized as legal (2864 mm FL) or sublegal size. Stepwise model selection based on Akaike information criterion 
(AIC) scores was used to evaluate the influence of these predictors. For each model, the AIC corrected for small 
sample sizes (AICc), the difference in AICc between the model and the best model (AAICc), and the likelihood ratio 
test statistic used to determine support of each model relative to a null model are provided. For each predictor in 
each model, the beta coefficient (8) and its standard error (SE) and the hazard ratio and its 95% confidence interval 
(95% CI) are provided. An asterisk (*) indicates that the predictor or model is significant (P<0.05). 


Hazard 


AlCce AAICe te ratio 


Full model 
44.88 11.28 
0.01 0.005 
—0.81 TES 
0.07 0.107 
-17.7 15,650 
1.09 1.506 
-16.7 86,490 
0.19 0.430 


Second-best model 
34.26 0.67 


Best model 
33.59 0.00 


Null model 
38.40 4.81 


Hazard Likelihood 
ratio 95% ratio Overall 
CI statistic ie 


1.00—1.02 
0.05-—4.47 
0.87—-1.32 
~0 to oo 
0.16—57.2 
~0 to oo 
0.52—2.82 


1.00-1.02 
0.43-85.1 


1.00—1.02 


Model with categorical predictor: fish length (sublegal or legal) 
30.75 


LG 3.01 1.120 20.22 


the first minute after release to the end of the acoustic 
monitoring period (Fig. 5). Bottom temperature, examined 
for internally tagged fish for which data were avail- 
able (n=35), was not a significant factor in PRM (Suppl. 
Table 7) (online only). Handling and fight times, examined 
for fish for which data were available (n=49), also were not 
significant factors in PRM (Suppl. Table 8) (online only). 


Discussion 


In this study, we found relatively high survivorship from 
recreational fishing methods (85%) for greater amberjack 
in the northern Gulf of Mexico. After accounting for all 
estimated sources of discard mortality, the probability of 
survival from PRM (89%) was lower than the probabil- 
ity of survival after AVM and CHM (95%). Notably, the 
recapture rate was relatively high (12%), underscoring 
the high fishing pressure in the area and likely increasing 


2.25-181 0.007* 


the chance of fish having repeated exposure to potential 
risks associated with discard mortality. In this study, fish 
size was the best predictor of PRM, and abiotic variables 
and factors associated with barotrauma (site depth) and 
mitigation of barotrauma (descender device use) were not 
associated with PRM. 

Results from Cox proportional hazards modeling with 
data for internally tagged fish (from 2019 and 2020) indi- 
cate a substantially (20 times) increased mortality risk 
for legal-size fish. The survivorship curve for legal-size 
and sublegal-size fish indicates that, from 7.6 h to 40 d 
postrelease, mortality probabilities were 2% for legal-size 
and 44% for sublegal-size fish. Therefore, when predicting 
postrelease outcome, it is important to consider that fish 
discarded by anglers during closed seasons may include 
large fish (=>864 mm FL), with a higher PRM risk than 
smaller fish, but during the open recreational fishing peri- 
ods anglers are more likely to discard sublegal-size fish 
when targeting legal-size fish for take. 
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Figure 5 


The Kaplan—Meier survivorship curve for greater amberjack (Seriola dumerili) 
of legal size (=864 mm in fork length [FL]) and sublegal size (<864 mm FL) 
with internal acoustic tags, (A) for the entire period (>50 d) for which transmit- 
ters and receivers were active and (B) for the first 38 h postrelease. Each cross 
along the curve indicates when data were censored because a fish emigrated 
away from acoustic receivers. The gray area represents the 95% confidence 
interval for survival probability. Fish were caught and tagged on charter boats 


in the northern Gulf of Mexico in 2019 and 2020. 


Having a body size smaller than legal size may reduce 
PRM. Results of exploratory Spearman rank correlation 
tests of predictor variables indicate a correlation of fight 
time with body size (Suppl. Table 1) (online only). This 
relationship is predicted because bigger fish are expected 
to fight harder. Although fight time was not explicitly 
predictive in Cox models, such a relationship may be con- 
founded by body size. Long fight times are expected to 
increase stress, elevate oxygen demand, and cause 
exhaustion, all of which could increase discard mortality 
risk (including risk of PRM). Fight and handling times 
have been reported to have physiological effects that 
reduce postrelease survival in blacktip sharks (Carchar- 
hinus limbatus) (Mohan et al., 2020). In our study, as 
expected, fish size also correlated negatively with han- 
dling time. These factors did not influence mortality in 
our Cox models but may be confounded by size. Effects of 
handling time may be less pronounced for fish discarded 
during normal recreational fishing compared to those in 
this study that involved tagging. Sublegal-size fish were 
more likely to be caught by jigging and to have observ- 
able fishing injuries, yet these factors did not appear to 
be correlated with PRM. Lastly, surface dissolved oxygen 
was a potential confounding variable in our study that 
weakly and negatively correlated with fish size in our 
data set, perhaps because more large fish were caught in 
summer in 2018 when most cases of low surface dissolved 
oxygen were observed in our study. 


10. 20930 
Time (h) 


Seasonally warm temperatures in 
autumn (~18°C) have been considered a 
likely factor of PRM in a cold-temperate 
species, the haddock (Melanogram- 
mus aeglefinus), in the Gulf of Maine 
(Capizzano et al., 2019). Cooler tempera- 
tures (e.g., at the surface, ~24°C in winter 
and spring versus 31°C in summer) have 
also been associated with increased post- 
release survival of red snapper in the Gulf 
of Mexico (Curtis et al., 2015). We did not 
find direct support for temperature as a 
contributor to PRM of greater amberjack 
in the Gulf of Mexico. Note, however, that 
bottom temperature data were limited 
in our study for fish sampled in spring 
in 2019, and the sample size of observed 
mortalities is small (n=7). Therefore, the 
differential between surface and bottom 
temperatures and thermocline depth 
should be considered as potential factors 
in future studies. 

One goal of our study was to assess 
the effects of 2 factors, site depth and 
descender device use, on PRM. We ini- 
tially aimed to attach acoustic teleme- 
try transmitters to external dart tags 
(Floy FH-69 tags). However, after a high 
rate of transmitters being shed in 2018 
(56% of fish shed their tags, confirmed 
by the reported recapture of 17% of 
tagged fish), we modified our approach to monitor PRM 
for a longer duration after fish release. Although tag shed- 
ding was frequent in 2018, data for detections of tags that 
were later shed indicate survival over the first several 
days after release: for the 9 fish inferred to have shed tags, 
tags appeared to remain attached between 2.8 and 33.2 d 
(median: 7.2 d), and 4 recaptured fish lived beyond the 
study period. Only 2 fish with external acoustic transmit- 
ters died after release, and these deaths occurred between 
2h and 2.8 d postrelease. Both of these fish were from rela- 
tively shallow sites (with depths of 33 and 37 m), with one 
released at the surface and one released with a descender 
device. Therefore, there was no evidence of barotrauma- 
induced PRM in 2018. In 2019 and 2020, we placed trans- 
mitters in the posterior end of the coelomic cavity and did 
not attempt to vent the swim bladder. Although no rup- 
tured swim bladders were observed during tag implan- 
tation, it is possible that some swim bladder gas was 
inadvertently released. For fish tagged in 2019 and 2020, 
there was no evidence that site depth increased mortality. 
Use of descender devices did not appear to mitigate PRM 
risk, but they did not add mortality risk, which might have 
been expected from predation of tethered fish on descent 
or from increased handling time when the descender 
device was used. 

Susceptibility to barotrauma varies among species and 
capture depths (Jarvis and Lowe, 2008). Species that 
ascend rapidly and undergo diel depth migrations are 
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predicted to resist depth-related trauma. In a recent study 
in which discard mortality in 2 species, the red snapper 
and gray triggerfish (Balistes capriscus), was examined, 
descender device use was found to be associated with a 
decrease in PRM for red snapper but not for gray trig- 
gerfish (Bohaboy et al., 2020). In addition, barotrauma 
risk increases with capture depth. In a recent study on 
PRM of greater amberjack at 2 sites that varied in depth, 
site depth was found not to predict survivorship (Jackson 
et al., 2018). In the Jackson et al. (2018) study, fish were 
vented prior to release, given that previous observations 
(Murie and Parkyn”) indicate turgidity in swim bladders 
of fish caught at depths greater than 45 m. Therefore, the 
potential for capture depth as a factor in postrelease sur- 
vivorship could still exist without mitigation of venting. 
Murie and Parkyn? noted that anatomy of greater amber- 
jack in the region of the pectoral girdle lacerates near the 
medial supracleithrum upon depth-related swim bladder 
expansion and allows bubble release. This “self-venting 
mechanism” (Murie and Parkyn’”), however, could still 
pose injury risk to released fish. 

In a mark-recapture study of greater amberjack that 
were collected at depths of 15—95 m off the Atlantic coast 
of the southeastern United States and were vented prior 
to release, no trends in recapture rates related to site 
depth were observed (i.e., no increased mortality occurred 
at sites where barotrauma would be predicted) (McGovern 
et al., 2005). In our study (in which fish were captured at 
depths <65 m), it was not known if fish mitigated baro- 
trauma by self-venting or if barotrauma simply was not 
a factor. It was also not known if greater amberjack faced 
higher PRM risk at capture sites with depths greater than 
64 m. Greater amberjack occuppied a wide depth range 
throughout the day. Even at deep reef sites (>64 m), fish 
may be hooked at shallower depths and require less com- 
pensation for luminal gas expansion of the swim blad- 
der than would be expected if they were hooked near the 
seafloor. 

Fish that succumbed to PRM in this study tended to die 
quickly (median: 1.2 h postrelease; range: 0.02—68.00 h 
postrelease). Sixteen acoustically tagged fish (20.5%) were 
not detected for more than 68 h (median: 0.2 h; range: 
0.00—3.50 h). The fate of these fish remains a question, 
although no evidence of predation existed. Emigration 
immediately following release is often reported from stud- 
ies in which passive acoustic telemetry was used, and as a 
result, the fates of many fish remain unknown, leaving the 
potential for higher PRM (Topping and Szedlmayer, 2011; 
Curtis et al., 2015). Two fish in our study were detected 
only once at the reef where they were tagged and released, 
and they were subsequently detected at a reef 4.5 km 
away. Therefore, initial emigration is clearly not always 


* Murie, D. J., and D. C. Parkyn. 2013. Preliminary release mor- 
tality of Gulf of Mexico greater amberjack from commercial and 
recreational hand-line fisheries: integration of fishing practices, 
environmental parameters, and fish physiological attributes. 
Southeast Data, Assessment, and Review SEDAR33-DW29, 
13 p. [Available from website.] 


associated with mortality. In our study, 16% of acoustically 
tagged fish that were detected alive >68 h postrelease 
moved between study sites with acoustic receivers, high- 
lighting that movement of greater amberjack among reefs 
is relatively common. Additional emigration in this study 
appeared associated with cyclonic storms. Six fish in 2020 
were not detected until 7-12 d after release. In addition, 
11 fish moved to different reefs that happened to have 
acoustic receivers, and 3 of those fish did not return to the 
reef where they were tagged and released. Therefore, hav- 
ing acoustic receivers in multiple locations helped prevent 
underestimation of postrelease survival and should be a 
consideration in studies of discard mortality in mobile reef 
fishes. 

Total discard mortality in our study, although compara- 
ble, was notably higher than the estimate for recreational 
fisheries (13.5% versus 10%, respectively) in the recent 
stock assessment of greater amberjack in the Gulf of Mexico 
(SEDAR, 2020). For greater amberjack in our study, PRM 
was the greatest component of discard mortality (AVM, 
CHM, and PRM). Discard mortality is an important compo- 
nent of total fishing mortality, and robust estimates of this 
rate are essential for effective stock management. 


Conclusions 


Results of this study, in which a broad range of abiotic fac- 
tors, such as site depth, use of a descender device, and fish 
length were examined, indicate that fish size presents the 
greatest risk for PRM of greater amberjack at sites with 
depths <64 m. Our study followed up previous work by 
Jackson et al. (2018) that examined discard mortality over 
a narrower depth range for fish that were vented prior 
to release. In our study, with fish released at depths of 
29-64 m, use of a descender device did not appear to influ- 
ence PRM. In contrast to what occurred in the Jackson 
et al. (2018) study, we found that PRM risk increases with 
fish size. These findings point to important management 
considerations, given that the release of small individu- 
als (<864 mm FL), required by restrictions on the size of 
fish that can be kept, appears to carry less PRM risk than 
the discard of legal-size fish. With our Cox proportional 
hazards model with fish size analyzed as a non-discrete 
variable, we estimated that mortality increases by 0.9% 
(hazard ratio: 1.009) per 1 mm increase in FL, a finding 
that should be considered when proposing changes to size 
regulations for greater amberjack in the northern Gulf of 
Mexico. 


Resumen 


EF] manejo de la pesca recreativa del medregal (Seriola 
dumerili) en el golfo de México incluye regulaciones de talla 
y temporadas de veda. La temperatura del agua, la salin- 
idad, los tiempos de forcejeo y manipulacion, asi como el 
barotrauma pueden influir en la supervivencia de los peces 
liberados. Examinamos la mortalidad y el comportamiento 
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posterior a la liberacion utilizando telemetria acustica de 
los movimientos de 78 peces, monitoreados hasta 58 dias a 
través de 3 esfuerzos de muestreo en 2018, 2019 y 2020 en 
un gradiente de profundidad de 29-64 m. Para evaluar los 
dispositivos de descenso como herramienta de mitigaci6én, 
asignamos peces a 2 tratamientos: liberaci6n en superficie 
sin ventilacién de la vejiga natatoria y liberaci6n con un 
dispositivo de descenso. Se utilizaron modelos de riesgos 
proporcionales de Cox para evaluar los efectos de la pro- 
fundidad del lugar, el tratamiento de liberacion, el tipo de 
carnada (sefuelo o carnada viva), las lesiones por pesca, 
las lesiones por marcado, los tiempos de forcejeo y manip- 
ulacién, las temperaturas de la superficie y del fondo, y 
la longitud de los peces. No se encontr6 una asociaci6n 
positiva entre el riesgo de mortalidad y la profundidad 
del sitio de liberacién, como se esperaria del barotrauma, 
ni una mayor supervivencia de los peces liberados con 
un dispositivo de descenso. El modelo mejor sustentado 
consideré tnicamente la longitud de los peces como fac- 
tor de mortalidad tras la liberacién; los peces de tamafio 
legal (=>864 mm de longitud furcal) tuvieron un riesgo de 
mortalidad 20 veces mayor que los peces mas pequenios. 
Nuestros resultados indican que, debido a las restricciones 
de tamanio, los peces de talla inferior a la legal liberados se 
enfrentan a un riesgo de mortalidad mucho menor que los 
peces de talla legal. 
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Abstract—The Dolly Varden (Salvelinus 
malma), an iteroparous, facultatively 
anadromous Pacific salmonid, displays 
diverse life history and migration pat- 
terns. Using otolith microchemistry, we 
inferred that individuals sampled in the 
Nushagak Commercial Salmon Fishery 
Management District, Alaska (sample 
size [n]=30; mean fork length=597 mm), 
had entered saltwater at ages 4-7 and 
were in their first (26%) or second 
(74%) summer at sea. Most (88%) of the 
fish that had spent 2 summers at sea 
migrated in consecutive years, but 2 of 
them skipped a migration to marine 
waters, remaining in fresh water for an 
additional year after migrating to sea in 
their first year. Only 15% of the individ- 
uals with 2 summers at sea had anadro- 
mous mothers. In contrast, conspecifics 
sampled on the other side of Bristol 
Bay in the Egegik Commercial Salmon 
Fishery Management District started 
migrating at an earlier age, migrated 
more often, and more often had anad- 
romous mothers. Together, these results 
highlight the differences in life history 
among Dolly Varden and indicate that 
freshwater rather than marine condi- 
tions influence life history patterns, 
at least for fish within the Nushagak 
District. 
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Salmonids exhibit a wide variety of 
migratory patterns, among and within 
species (Klemetsen et al., 2003; Quinn, 
2018). Some species reside in small 
streams, whereas others migrate 
between streams and larger rivers 
(Northcote, 1997), migrate between 
streams and lakes (Arostegui and Quinn, 
2019), or are anadromous and spend 
weeks, months, or years at sea (Quinn 
and Myers, 2004). Anadromous salmo- 
nids that are iteroparous may vary in 
age of first marine migration (Swanson 
et al., 2010; Bond et al., 2015) and in the 
number or sequencing of those migra- 
tions (Spares et al., 2015). This variation 
is especially marked in the char species 
(Salvelinus spp.) (Dunham et al., 2008), 
with anadromy positively associated with 
key life history traits, such as growth 
rate, age at first spawning, fecundity, and 
longevity (Tallman et al., 1996), but also 
with high mortality (Brown et al., 2019). 

The Dolly Varden (S. malma) is an 
iteroparous char species of the North 
Pacific Rim with 8 subspecies (Taylor 
and May-McNally, 2015), and there are 


fluvial (Denton et al., 2010; Ayer et al., 
2018), adfluvial (Markevich et al., 
2018), and anadromous or partially 
anadromous populations of this species 
(Bond and Quinn, 2013; Morrison et al., 
2021). Anadromy is known for the 2 
North American subspecies or forms 
in Alaska. The southern form of Dolly 
Varden, S. m. lordii, found south of the 
Alaska Peninsula, inhabits estuarine 
and nearshore habitats (Bond et al., 
2014a) and may overwinter in marine 
waters (Bernard et al., 1995). In con- 
trast, the northern form, S. m. malma, 
found north of the Alaska Peninsula, 
makes more extensive marine migra- 
tions (Morita et al., 2009; Courtney 
et al., 2018) and typically overwinters 
in fresh water after spawning and then 
migrates to sea annually; however, 
individuals of this subspecies may skip 
a year between migrations, remain- 
ing in freshwater habitats (Gallagher 
et al., 2018). 

In Bristol Bay, in southwestern 
Alaska, populations of the northern form 
of Dolly Varden may be non-anadromous 
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(Denton et al., 2010) or partially anadromous, as 
has been inferred from catches at weirs (Lisac’) 
and from microchemical analysis of otoliths taken 
from fish sampled in marine waters at the mouth 
of the Egegik River (Hart et al., 2015) and taken in 
fresh water (Scanlon, 2000). Herein, we report life 
history data (size, total age, age at first migration, 
migration frequency, and maternal anadromy) on 
anadromous individuals of the northern subspe- 
cies of Dolly Varden sampled from the Nushagak 
Commercial Salmon Fishery Management District 
(hereafter, referred to as the Nushagak District), 
a commercial fishing zone within the Nushagak 
Bay estuary that opens into Bristol Bay. We com- 
pare them to data on conspecifics similarly sam- 
pled from the Egegik Commercial Salmon Fishery 
Management District (hereafter, referred to as the 
Egegik District) in the northeastern arm of Bristol 
Bay (Hart et al., 2015). The fish have access to 
common feeding areas, with the mouths of the riv- 
ers that drain into these bays only ~70 km apart 
(Fig. 1), allowing us to assess differences that might 
result from their respective freshwater ecosystems 
in the same region. 


Materials and methods 


On 25 June 2018, 30 fish identified as the north- 
ern form of Dolly Varden (Taylor and May- 
McNally, 2015) by. E. Taylor”, of the University 
of British Columbia, were taken as bycatch in 
the commercial gillnet fishery for sockeye salmon 
(Oncorhynchus nerka) in the Nushagak District 
(Salomone et al.°). The use of the samples exam- 
ined in this study was permitted by the Alaska 
Department of Fish and Game. These fish could 
not be sold as food and would have been dis- 
carded if they had not been made available to us 
for study. The fish were frozen and later thawed, 
measured (fork length [FL]), and weighed. Paired 
t-tests were used to detect statistically significant 
differences between fish from the Nushagak and 
Egegik Districts. 

Otoliths were removed for determination of life history 
characteristics by using visual and microchemical analyses. 
These analyses followed procedures used on conspecifics 
also taken as bycatch in the gillnet fishery between 22 June 
and 2 July 2013 in the Egegik District, on the opposite side 
of Bristol Bay from the Nushagak District (Hart et al., 


' Lisac, M. J. 2010. Abundance and run timing of Dolly Varden 
in the Middle Fork Goodnews River, 2008 and 2009. Alsk. Fish. 
Data Ser. Rep. 2010-13, 14 p. [Available from Togiak Natl. Wildl. 
Refug. Off., 6 Main St., Dillingham, AK 99576-0270.] 

2 Taylor, E. 2018. Personal commun. Dep. Zool., Univ. British Colum- 
bia, Beaty Biodivers. Cent., Rm. 310, Vancouver, Canada V6T 124. 

3 Salomone, P., T. Elison, T. Sands, J. Head, and T. Lemons. 2019. 
2018 Bristol Bay annual management report. Alsk. Dep. Fish 
Game, Fish. Manag. Rep. 19-12, 64 p. [Available from website. ] 
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Figure 1 


Map of the 2 locations in southwestern Alaska where individuals of 
the northern subspecies of Dolly Varden (Salvelinus malma malma) 
were taken as bycatch in the commercial gillnet fishery for sock- 
eye salmon (Oncorhynchus nerka). Samples were collected from the 
Nushagak Commercial Salmon Fishery Management District on 
25 June 2018 for this study, and comparable samples of the same 
subspecies of Dolly Varden were collected in the Egegik Commercial 
Salmon Fishery Management District by Hart et al. (2015) on 22 
June and 2 July 2013. 


2015). The similarities in gillnet mesh size (typically a 
mean diameter of 12.7 cm [standard deviation (SD) 0.3]) 
and in the times when these fisheries are open made it pos- 
sible for us to compare life history patterns between these 
2 collections. It is important to note, however, that neither 
collection fully represents the population in the district 
from which it was sampled owing to the limited collection 
dates and large basins from which these specimens may 
have originated, basins that include not only the basins 
of the Nushagak and Egegik Rivers but also other basins 
throughout Bristol Bay and southwestern Alaska. 

Otoliths were extracted, sectioned, mounted, polished to 
expose the core and surrounding matrix, and aged visually 
by using reflected light. A subset of otoliths (sample size 
[n]=23) was selected for microchemical analysis on the 
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basis of the clarity of the preparation and ease of deter- 
mination of the otolith’s core from the surrounding area. 
Strontium (Sr) and barium (Ba) are incorporated into fish 
bony structures relative to their concentrations in the 
surrounding water. The concentration of these elements 
is typically correlated with salinity (Elsdon et al., 2008) 
and is highly conserved in fish otoliths (Kraus and Secor, 
2004). Normalized to calcium (Ca) levels, Sr and Ba con- 
centrations in otoliths have been used to indicate migra- 
tory patterns and life history in diadromous fish species 
(Walther, 2019) and especially in salmonids (Swanson 
et al., 2010; Bond et al., 2015; Austin et al., 2019). 

Laser-ablation mass spectroscopy of otoliths was con- 
ducted with an iCap RQ ICP-MS‘ (Thermo Fisher Scien- 
tific, Waltham, MA). Element concentrations, adjusted to 
NIST-612 standards, were converted to molar mass ratios 
of *Sr:*°Ca and '°°Ba:*°Ca, presented as millimoles per 
mole and micromoles per mole, respectively. Winter annuli 
were visually identified by using dark growth rings and 
were measured from the core along the ablated transect. 
These measurements were then converted to a ratio of 
the total transect measured, and the corresponding ratio 
of lengths along the transect were overlaid onto a 9-point 
rolling mean of Sr:Ca and Ba:Ca ratios for individual oto- 
liths to infer initial age of migration and annual migra- 
tory and overwintering patterns. Maternal anadromy was 
identified when core Sr:Ca ratios were elevated above the 
surrounding freshwater signals, reflecting the concentra- 
tions in the female during egg production (Kalish, 1990). 
Similarly, anadromy was inferred when the Sr:Ca ratio 
was >3 mmol/mol and the Ba:Ca ratio was <6 pmol/mol 
(Bond et al., 2015). In all cases, there was agreement 
between Sr:Ca and Ba:Ca ratios, indicating anadromy. 
The life history and microchemical data from our study 
can be accessed from GitHub. 


Results and discussion 


The Dolly Varden sampled from Nushagak District had 
an average length of 597.1 mm FL (SD 40.6), with a 
range of 522-690 mm FL, and an average mass of 2.16 kg 
(SD 0.42), with a range of 1.63—3.25 kg. Both of those 
means are greater than those of the fish caught in the 
Egegik District at essentially the same time of year, 
with a mean length of 509.0 mm FL (t=10.0, P<0.001) 
and a mean mass of 1.49 (t=4.1, P<0.001). However, the 
mean age of 6.1 years (range: 5-8 years) for fish from 
the Nushagak District did not differ from the mean age 
of 5.7 years (range: 4—9 years) for fish from the Egegik 
District (t=1.75, P=0.09) (Suppl. Table) (online only). These 
length and age values are comparable to those for con- 
specifics at similar latitudes (Maekawa and Nakano, 
2002), although older individuals occur in Arctic popula- 
tions (Stolarski and Sutton, 2013; Gallagher et al., 2016). 


4 Mention of trade names or commercial companies is for identi- 
fication purposes only and does not imply endorsement by the 
National Marine Fisheries Service, NOAA. 


Morrison et al. (2021) reported that rapid growth in 
fresh water was associated with early age at first migra- 
tion in Dolly Varden and that growth rate in fresh water 
is typically inversely related to smolt age in salmonids 
(McCormick et al., 1998; Quinn, 2018). In contrast to the 
modest differences in size and age between the popula- 
tions on opposite sides of Bristol Bay, a clear difference in 
age at migration was observed. The fish sampled from the 
Nushagak District first migrated to sea much later in life 
than the fish sampled from the Egegik District, with mean 
ages at migration of 5.2 years (range: 4-7) and 2.5 years 
(range: 2-3), respectively (t=15.33, P<0.001), indicating 
that the fish from the Nushagak District grew at a slower 
rate than the fish from the Egegik District. Similarly, a 
population of the southern subspecies of Dolly Varden 
along the southern Alaska Peninsula also migrated to 
sea at a younger age, although at a smaller body size at 
that age of migration (Bond et al., 2015). Dolly Varden at 
the northern end of their range also smolt at a younger 
age than the fish sampled from the Nushagak District 
(Morrison et al., 2021). 

In addition to the difference between populations in age 
at first migration, 6 of the fish from the Nushagak District 
were in their first summer at sea, and the other 17 fish were 
in their second summer (mean age: 1.7 years) compared 
with the mean of 4.2 summers at sea (range: 3-8 summers) 
for the fish from the Egegik District. 

Skipped migrations (i.e., remaining in fresh water 
for 2 winters rather than 1 winter between migrations) 
(Fig. 2A) were not common. Examining fish with multi- 
ple seasons at sea, we found that only 2 of 17 fish from 
the Nushagak District had skipped a migration to sea 
and that no fish from the Egegik District had skipped a 
marine migration. Skipped migrations were much more 
common (76% of the sample) in a study of fish sampled 
from an Arctic river that had a long (~800 km roundtrip) 
migration route (Gallagher et al., 2018). In comparison, 
the mouth of the Egegik River is only 45 km from Lake 
Becharof. Individuals from the Nushagak District migrate 
an unknown distance but likely less than the 450-km 
length of the Nushagak River because the Dolly Varden is 
common in the basin closer to Bristol Bay (May-McNally 
et al., 2015). 

Finally, the incidence of maternal anadromy in our 
sampled individuals differed greatly (e.g., as evident in 
the data for 2 specimens shown in Figure 2, A and B): 
only 3 of 20 fish from the Nushagak District had mothers 
that migrated to sea in the season before spawning com- 
pared with 14 of 25 fish from the Egegik District having 
such mothers, yet all samples were collected in marine 
waters. 

A skipped migration would create a freshwater chemi- 
cal signal in a female that would not be distinguishable 
from a non-anadromous fish. In our sample, skipped 
migrations were rare; therefore, our results indicate that, 
in both populations, maternal migrations may not be com- 
mon. However, a sample of juvenile Dolly Varden from 
tributaries of Lake Aleknagik, in the Wood River system 
of the Nushagak River basin, contained only individuals 
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Figure 2 
Profiles of the ratio of strontium to calcium (Sr:Ca) for otoliths from 2 Dolly 
Varden (Salvelinus malma) sampled on 25 June 2018 from the Nushagak 
Commercial Salmon Fishery Management District in southwestern Alaska, 
truncated to show the core at the origin. The vertical gray solid lines indicate 
the approximate positions of annuli. The horizontal dark gray dashed lines 
indicate the Sr:Ca threshold for interred marine residence (3 mmol/mol). 
(A) The top profile indicates maternal anadromy (elevated Sr:Ca at the core) 
and a year in freshwater between anadromous migrations (i.e., a skipped 
migration) at ages 4+ and 6+. (B) The bottom profile has no indication of 
maternal anadromy and indicates capture during the initial migration of the 
fish at age 7+. For clarity of the graphs, values of the ratio of barium and Ca 


been more common in the absence of 
fishing. We also acknowledge that inter- 
preting age and microchemical data from 
otoliths is not without error (Stolarski 
and Sutton, 2013). For example, expos- 
ing otolith cores can be less reliable for 
adults than for juveniles, diluting or 
masking the core signature (Volk et al., 
2000). 

In addition, we made conclusions with 
the assumption that fish originated 
from the natal rivers within the district 
where they were caught. Dolly Varden 
have strong homing instincts, as indicated 
by results from tagging (Armstrong, 1974) 
and population genetics (Bond et al., 
2014b) studies, but some individuals, 
especially those of the northern subspe- 
cies, can disperse widely (DeCicco, 1992; 
Morita et al., 2009). Still, the differences 
in life history observed in the samples 
from the Egegik and Nushagak Districts 
indicate that the examined samples did 
not represent a single mixed population 
but, rather, ones that were segregated 
to at least some extent. Despite issues 
with methods, common to the sampling 
of fish in both the Nushagak and Ege- 
gik Districts, and despite the access that 
these fish have to common, proximate 
foraging areas, the dramatic differences 
in life history of these populations indi- 


are not shown. 


whose mothers had been to sea in the summer prior to 
spawning (Dennert et al., 2016). The apparent differences 
among these results likely reflect the range of life history 
patterns of Dolly Varden in these large basins and the fact 
that our samples likely did not fully capture this varia- 
tion. These variable life history patterns likely arise from 
differential survival and reproductive success, and there- 
fore from differences in life history traits, across the het- 
erogeneous freshwater landscape within the Nushagak 
River basin. These varied freshwater conditions may have 
more influence on the survival and reproductive success of 
fish migrating from the basin than the marine conditions 
to which all fish that migrate to sea are exposed. Regard- 
less of the causes, however, the results indicate that the 
life history of the mother often differs from that of her 
offspring. 

The fish used in this study were taken in gill nets that 
are size-selective (Todd and Larkin, 1971), affecting the 
range of phenotypes observed in a sample. Moreover, fish- 
ing mortality in general reduces age and longevity in pop- 
ulations (Sharpe and Hendry, 2009). Therefore, we cannot 
infer that fish smaller than the fish in our study do not 
enter marine waters, and repeat migrations might have 


cate a dominant role of freshwater fac- 
tors rather than marine factors on these 
traits. In this regard, the life history 
differences echo those found among pop- 
ulations of other salmonid species, including the sockeye 
salmon of Bristol Bay (Quinn et al., 2009). 


Resumen 


El salmon Salvelinus malma, es una especie iter6para 
del Pacifico, facultativamente anddroma que presenta 
ciclo vital y patrones de migracion diversos. A partir de la 
microquimica de los otolitos, deducimos que los individuos 
muestreados en el Distrito de Administracion de la pes- 
queria comercial de salmén de Nushagak, Alaska (tamano 
de muestra [n]=30; longitud furcal promedio=597 mm), 
habian entrado en el agua salada a los 4—7 afios y estaban 
en su primer (26%) o segundo (74%) verano en el mar. La 
mayoria (88%) de los peces que habian pasado 2 veranos 
en el mar migraron en afios consecutivos, pero 2 de ellos 
se saltaron una migraciOn a aguas marinas, permaneci- 
endo en agua dulce un afio mas después de haber migrado 
al mar en su primer ano. Sdlo el 15% de los individuos 
con 2 veranos en el mar tenian madres andédromas. Por 
el contrario, los congéneres muestreados al otro lado de 
la bahia de Bristol, en el Distrito de Administracién de 
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la pesqueria comercial de salmon de Egegik, empezaron 
a migrar a menor edad, con mas frecuencia, y tuvieron 
mas a menudo madres anddromas. En conjunto, todos 
estos resultados ponen de manifiesto las diferencias en el 
ciclo vital de este salmon e indican que son las condiciones 
de agua dulce y no las marinas las que influyen en los 
patrones del ciclo vital, al menos en el caso de los peces del 
distrito de Nushagak. 
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Abstract—Crustacean species are 
socioeconomically and ecologically cru- 
cial across the world. For crustaceans, 
as ectotherms, anthropogenic climate 
change threatens to significantly alter 
key life history characteristics, such 
as size at maturity and growth rate. 
Because crustaceans are difficult to 
age, length data are used in assess- 
ments of crustacean stocks; however, 
climate-induced changes in maturation 
and growth can greatly influence the 
performance of size-structured stock 
assessment models. We _ coupled 
individual-based and _ size-structured 
models for American lobster (Homarus 
americanus) off northeastern North 
America in the Gulf of Maine—to 
conduct a novel sensitivity analysis 
of the effects of maturity and growth- 
related input parameters on model 
outputs. For this analysis, we used 
a bottom-up approach (with param- 
eters shifted independently) and a 
top-down approach (with parameters 
shifted jointly as they were predicted 
to be influenced by climate change). 
We found that our American lobster 
stock assessment model is resilient to 
relatively extreme shifts in biological 
input parameters. For size-structured 
modeling in assessments of crustacean 
stocks, we recommend the expansion of 
sensitivity analyses to include evalua- 
tion of the influence of climate-driven 
changes on input parameters based on 
time-varying life history traits. 
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Anthropogenic climate change is trans- 
forming many marine’ ecosystems 
through warming waters, ocean acidi- 
fication, freshening, and deoxygenation 
(Brander, 2010; Doney et al., 2012; 
Gattuso et al., 2018; Doney et al., 2020). 
Perturbations to the abiotic environ- 
ment, in particular to temperature, 
are especially influential on marine 
ectotherms because they do not physi- 
ologically regulate their body tempera- 
ture (Huey et al., 1993; Madeira et al., 
2012); rather, their body temperature 
is driven by the environment. As a con- 
sequence, temperature directly influ- 
ences individual and population-level 
biological processes of crustacean spe- 
cies, such as metabolism, recruitment, 
reproduction, growth, size at maturity, 
and natural mortality (Madeira et al., 
2012), which have significant implica- 
tions for assessment and management 
of crustacean fisheries (Audzijonyte 
et al., 2016). Typically, in data-rich 
crustacean stock assessments, size- 
structured models are used (Punt et al., 
2013), and the outputs of such models 
can be influenced by environmentally 
driven variability in input parameters 
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Fisheries and Oceans Canada 
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Nanaimo, British Columbia V9T6N7, Canada 


based on size-related life history char- 
acteristics, such as growth and size 
at maturity. Therefore, it is important 
to quantify how climate-driven shifts 
in key life history characteristics will 
influence crustacean stocks and will 
manifest in assessment models used for 
guiding future management decisions. 

The American lobster (Homarus 
americanus) is an ecologically and socio- 
economically vital crustacean species in 
the northwestern Atlantic Ocean (Le 
Bris et al., 2018), and the biology of this 
species is directly influenced by tem- 
perature. American lobster, like many 
crustaceans, grow through a series of 
molts, also known as ecdysis. During 
ecdysis, the old carapace is replaced 
with a new, larger one (Herrick, 1911). 
Molting typically occurs annually in 
adult American lobster, although it can 
happen more than once or be skipped 
entirely depending on the size, age, and 
maturity of the individual (Herrick, 
1911; Aiken and Waddy, 1976; Aiken, 
1977; Comeau and Savoie, 2001). 
Although the physiology of individual 
American lobster is known to influ- 
ence growth processes, temperature is 
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a primary abiotic driver of growth changes in this species. 
Rising temperatures have been shown to increase molting 
frequency and decrease molting increment, the length a 
lobster grows in a given molting event (Aiken, 1977). Addi- 
tionally, in several studies, rising temperatures were found 
to contribute to reduced sizes at maturities for American 
lobster (Little and Watson, 2003, 2005; Le Bris et al., 2017; 
Waller et al., 2021). Indeed, climate-driven changes in these 
life history characteristics can likely affect results from the 
size-structured stock assessment model currently used for 
management of American lobster (ASMFC, 2020). 

Understanding the effect that climate-driven shifts in 
input parameters for life history characteristics have on 
stock assessment model outputs is critical for future model 
development and fisheries management. In a recent study, 
incorporating temperature-driven recruitment was found 
to improve the performance of a size-structured stock 
assessment model for American lobster (Tanaka et al., 
2019). When simulating the effects of pooling multiple pop- 
ulations of red rock lobster (Jasuws edwardsii) with varying 
growth rates, the performance of a size-structured stock 
assessment model was not reduced (Punt, 2003), indicat- 
ing that accounting for different growth rates of assessed 
populations may not be consequential for estimating ref- 
erence points. If pooling population data of red rock lobster 
had reduced model performance, it may have indicated a 
need to further consider the importance of variable growth 
in future assessments of the population. In contrast, other 
research results indicate that failing to account for the 
plasticity of growth in fisheries stock assessment models 
can lead to deviations of more than 30% in outputs, criti- 
cally altering the calculation of reference points (Lorenzen, 
2016). Indeed, depending on the species biology and stock 
assessment model design, changes in growth can have 
inconsistent effects on model outputs. 

Typically, sensitivity analyses can evaluate whether 
uncertainties in model assumptions, input data sources, 
and biological parameters have an effect on reference points 
or other model outputs (Maunder and Punt, 2013; Maunder 
and Piner, 2015). However, in these analyses, adjustments 
to inputs usually are considered only on their own, rather 
than in combination (Lehuta et al., 2010; Saltelli et al., 
2019), and these analyses seldom involve testing model sen- 
sitivity to inputs that are based on life history traits and 
developed outside of the assessment model, such as growth 
transition matrices. Given the potential for dissimilar con- 
sequences of changing life history characteristics for stock 
assessment model outputs and the yet unrealized but possi- 
ble shifts in crustacean growth in the future, it is important 
to evaluate the sensitivity of size-structured stock assess- 
ment models on a case-by-case basis. 

In our study, we conducted a novel sensitivity analysis 
of a length-structured stock assessment model for Ameri- 
can lobster by using an individual-based simulation model 
to evaluate the sensitivity of the model to shifts in growth- 
related life history input parameters, specifically molting 
probability, molt increment probability, and size at matu- 
rity. In a series of sensitivity analyses, we used classical 
bottom-up methods in which each parameter was shifted 


independently, and we used a top-down approach in which 
parameters were jointly shifted under the driving mecha- 
nism of climate change. We used both approaches to ascer- 
tain at what point shifts in these input parameters could 
result in a significant change in the stock status time 
series estimated with the length-structured stock assess- 
ment model, relative to a historical baseline. Our overar- 
ching goal for this study was to determine the extent of the 
influence that climate change has on the reliability and 
robustness of outputs from the stock assessment model 
used for American lobster in the Gulf of Maine. 


Materials and methods 
Shifting growth and size at maturity 


Seasonal growth matrices in this study were calculated by 
using an individual-based lobster simulator model (IBLS) 
first developed by Chen et al. (2005) and later expanded 
by Chang (2015) and Mazur et al. (2018). This model sim- 
ulates individual lobster from recruitment to mortality 
by sending each lobster through random Bernoulli tri- 
als representative of life history and fishery parameters 
derived from prior field research and modeling (Chen 
et al., 2005; Chang, 2015; Mazur et al., 2018). This sea- 
sonal probabilistic model is used to simulate lobster fish- 
ery dynamics to capture complex fishery-dependent and 
fishery-independent processes (Chen et al., 2005; Zhang 
et al., 2011) and has historically been used to test the per- 
formance of the American lobster stock assessment model 
(Chen et al., 2005). The model creates individual lobster 
records over a given time series that include sex, size bin, 
carapace length (CL), maturity, and mortality, allowing 
calculation of population abundance, spawning stock bio- 
mass, and landings (Mazur et al., 2018). A full explanation 
of this model can be found in Mazur et al. (2018). 

The IBLS can be used to create seasonal growth matri- 
ces by simulating lobster with total absence of fishery- 
dependent and fishery-independent mortality as well 
as recruitment. Effectively, the abundance of lobster 
remains constant over the simulated time series, but the 
biomass changes exclusively because of input of data on 
growth of the animals. At each step, a lobster is in 1 of 35 size 
bins (in increments of 5 mm, from 53 mm CL to 2223 mm 
CL). The simulation was run long enough that every lobster 
ended up in the final size bin at the end of the time series. 
Every growth instance for every lobster for a given season 
over the entire time series is marked in a matrix with the 
size bin before the molting event on the x-axis and the size 
bin after the molting event on the y-axis, and the matrix is 
scaled so that the sum of each row is effectively 1. This sim- 
ulation creates a probabilistic growth matrix in which each 
row is a function of size change for a given lobster of that 
size class. This process is done 4 times, once for each sea- 
son: winter (January—March), spring (April—June), summer 
(July-September), and fall (October—December). 

The growth inputs to the IBLS are 2 independent fac- 
tors: molting probability and probability for different molt 


242 


Fishery Bulletin 120(3—4) 


increments. Molting probability is the probability of a lob- 
ster molting in a particular time step dependent on the 
CL, maturation status of the individual, and how many 
seasons it has gone without molting (Fig. 1). Molt incre- 
ment probability is the probability of a lobster growing a 
certain size (1-20 mm CL in 1-mm-CL bins) because of a 
molting event and is dependent on the CL of the individ- 
ual (Fig. 1). The input data for these parameters in the 
base case of the IBLS came from the stock assessment con- 
ducted in 2015 (ASMFC, 2015). 

Because of climate change, American lobster are 
expected to molt more frequently but grow less per molt 
(ASMFC, 2015). To simulate these effects on overall 
growth, both molting probability (P),) and molt increment 
probability (Py) were manipulated in the IBLS. Molting 
probability was increased by shifting left in relation to 
seasons since the prior molt (Fig. 1) and was described by 
the following equations: 


as + 
Py = vor, (1) 
Rez, 
1,...,k¢j, ifimmature 
yas = (2) 
2,...,kcy, if mature; and 
pope aes o~8-08127 + 0.076535CL (3) 


where yas = time spent (in units of the time step of the 
model; in this case, the unit is seasons) at cur- 
rent size of an individual lobster; 
kcoy = the longest time a lobster of a given CL 
(in millimeters) could feasibly go before molt- 
ing (NEFSC, 1996; ASMFC’); and 
6, = the shifting parameter. 


Therefore, a 6, of 1 would represent a shift of 1 season, 
increasing the overall probability of molting in compari- 
son to the unshifted probability. 

Average size increase per molt was lowered by shift- 
ing molt increment probability left in relation to the size 
increase per molt (Fig. 1), described by using the equations 
below: 


Pyn = N(ALy, — by,0”), (4) 


1.2236 +0.1294L for males with L < 95; 
— /|1.2236+ (0.1294 x95) for males with L = 95; (5) 


1.2288 + 0.1285L for females with L < 82; and 
1.2288 + (0.1285 x 82) for females with L > 82, 


where N =the normal distribution truncated by upper 
and lower boundary probabilities of 0.975 and 


1 ASMFC (Atlantic States Marine Fisheries Commission). 2000. 
American lobster stock assessment report for peer review. Atl. 
States Mar. Fish. Comm., Stock Assess. Rep. 00-01 (Suppl.), 
334 p. ASMFC, Arlington, VA. [Available from website.] 


0.025, respectively, with standard deviation (0) 
being equal to 2.1 (ASMFC”); 
__L =the current length (in millimeters); 


AL,, = the change in length (in millimeters), given L 
and sex; and 
b. = the shifting parameter. 


Therefore, a b, of 1 would represent a shift of 1 mm, 
decreasing the overall size increment change during a 
given molt. 

To maintain some biological realism, shifts of molting 
probability and molt increment probability were paired, 
and the corresponding growth matrices reflect possible 
effects from climate change. Two paired shifts were used 
in this study. The first paired shift, hereafter referred to as 
G1, was a leftward shift of molting probability by 1 year 
and of molt increment probability by 1 size bin (6,=6,=1) 
(Fig. 1). The second paired shift, hereafter referred to as 
G2, was a leftward shift of molting probability by 2 sea- 
sons and of molt increment probability by 2 size bins 
(6,=b5=2) (Fig. 1). 

Probabilistic size at maturity (Pg) in the IBLS is cal- 
culated with the below equation: 


1 
Fsam = 14 003% (CLL) ’ (6) 
where Pgay = the probability of maturity of an individual 
lobster of a given CL (in millimeters); and 
Lz, = the predefined CL (in millimeters) at 50% 
maturity. 


The parameter L;,. was set to 90.81 mm CL for the IBLS 
base case (ASMFC, 2015). Given that size at maturity for 
American lobster is expected to decrease 2.8 mm CL per 
1°C rise in bottom temperature (Le Bris et al., 2017) and 
that current projections of bottom temperature for the 
Gulf of Maine are for an increase up to 2°C by 2050 and up 
to 4°C by 2100 (IPCC, 2019; Brickman et al., 2021), the L;5 
values of 85.21 and 79.61 mm CL were additionally tested 
in this study. These lengths are not considered projected 
values of size at maturity, but rather they are the L;, val- 
ues needed in calculations to provide reasonable levels of 
change from historical data for the sensitivity analyses. 
The IBLS was used to generate a total of 7 sets of growth 
matrices in this study (with 4 matrices in each set corre- 
sponding to seasons) (Table 1). The first set (the IBLS base 
case) was calculated with the original (unshifted) molt 
probability and molt increment probability paired with the 
original L;, value of 90.81 mm. The next 4 sets of matrices 
were calculated under IBLS scenarios 2—5, with shifts of 
either growth (G1 or G2) or Ls, (85.21 or 79.61 mm CL), 
and the last 2 sets were calculated under IBLS scenarios 6 
and 7, with paired shifts of both growth and Ls, (G1 and 
85.21 mm CL or G2 and 79.21 mm CL). Scenarios 2—5 were 


2 ASMFC (Atlantic States Marine Fisheries Commission). 2006. 
American lobster stock assessment report for peer review. Atl. 
States Mar. Fish. Comm., Stock Assess. Rep. 06-03 (Suppl.), 
175 p. ASMFC, Boston, MA. [Available from website.] 
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Figure 1 


The relationships (A) of molting probability in the summer to number of previous seasons with no molting 
and (B) of cumulative molt increment probability to size increase for an American lobster (Homarus ameri- 
canus) with a carapace length (CL) of 130 mm. The lines represent the probabilities in the individual-based 
lobster simulator model (IBLS) base case (or original scenario) and in each of the paired shifts from the 
IBLS base case used to simulate the effects of climate change on overall growth, a leftward shift of molting 
probability by 1 year and molt increment probability by 1 size bin of 1 cm CL (G1) and a leftward shift of 


molting probability by 2 seasons and molt increment probability by 2 size bins (G2). 


Table 1 


Individual-based lobster simulator (IBLS) scenarios for American lobster (Homarus ameri- 
canus) in the Gulf of Maine. The parameters of the scenarios are a paired set of shifts in molt- 
ing probability and molt increment probability for the growth dynamic and the predefined 
carapace length at 50% maturity (L;,) for each scenario. The IBLS base case or original sce- 
nario includes input data from the American lobster benchmark stock assessment conducted 
by the Atlantic States Marine Fisheries Commission in 2015. The paired shifts in growth 
factors are a leftward shift of molting probability by 1 year and molt increment probability 
by 1 size bin (G1) and a leftward shift of molting probability by 2 seasons and molt increment 


probability by 2 size bins (G2). 


Parameter 1 2 


Growth dynamic Original Gl 
Ls. (mm) 90.81 90.81 


used to observe effects from specific parameters, whereas 
scenarios 6—7 were meant to provide outputs that are more 
biologically realistic and expected given the predicted rela- 
tionships between climate change and the life history traits 
on which these parameters are based. Additionally, an 
analysis of time-varying growth and _ size-at-maturity 
parameters was conducted. Results of this analysis can be 
found in Supplementary Figure 1 (online only). 


Stock assessment and sensitivity analyses 


The University of Maine Lobster Stock Assessment Model 
(UMM) was initially developed by Chen et al. (2005) and 


IBLS scenario 


4 5 6 7 
Original Original Gl G2 
85.21 79.61 85.21 79.61 


expanded in ASMFC (2015) and Tanaka et al. (2019). It 
is a seasonal, sex-specific, length-structured assessment 
model for American lobster in the Gulf of Maine, Georges 
Bank, and Southern New England. It was designed with 
input from the Atlantic States Marine Fisheries Commis- 
sion with the intent that it would be used for future lob- 
ster stock assessments. The population dynamics equation 
of the UMM is as follows: 


PY M, 


Nig =|Nis-1 X Gs Xe + Ry; (7) 


where WN, , = a vector of the number of lobster in each size 
bin in year ¢ and season s; 
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G = the seasonal growth transition matrix; 

F = the seasonal fishing mortality; 

M = the seasonal natural mortality; and 

Re=recruitment abundance for each size bin 
(Chen et al., 2005). 


A list of all data used in the UMM consistently across sce- 
narios can be found in Table 2. For more detailed descrip- 
tions of this model, see Chen et al. (2005) and Tanaka et al. 
(2019). 

In the base case of the UMM, the original growth matri- 
ces and size at maturity of 90.81 mm CL from the IBLS base 
case were used as input data. Growth transition matrices 
and size-at-maturity data from the other 6 IBLS scenarios 
were individually input into the UMM, for a total of 7 sce- 
narios including the UMM base case (all growth matrices 
used are provided in the Supplementary Material [online 
only]). For each scenario, biological reference points (BRPs) 
were calculated for output reference abundance (individu- 
als larger than 53 mm CL) by using the methods outlined 
in ASMFC (2015): the target was calculated as the 75th 
percentile of reference abundance over the time series, and 
the threshold was calculated as the 25th percentile of ref- 
erence abundance over the time series (ASMFC, 2015). It is 
important to note that the reference time series for these 
calculations used by the Atlantic States Marine Fisheries 
Commission was for the period of 1982-2003, but in this 
study, we used data from a reference period of 1984—2003 
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because of data input limitations. These BRPs allowed 
determination of historical fishery status over time, mean- 
ing simply the reference abundance of a given year in rela- 
tion to the predefined BRPs (below the 25th percentile, 
between the 25th and 75th percentile, or above the 75th 
percentile). By using these reference points, terminal-year 
stock status was compared between all UMM runs in this 
study. However, for all sensitivity analyses in this study, 
historical fishery statuses over the entire time series were 
compared between each UMM scenario and the UMM 
base case. 

Scenarios 6 and 7 in the IBLS were designed to repre- 
sent small and large future climate effects, respectively. 
These effects on growth and size at maturity are plausi- 
ble given climate projections (IPCC, 2019; Brickman et al., 
2021), but it is unknown if these changes are large enough 
to affect stock status from what it would be under the 
UMM base case. To this end, a sensitivity analysis was 
conducted for scenarios 6 and 7 to examine whether stock 
status differed in consecutive years between the scenario 
and the UMM base case. 

In the sensitivity analyses, IBLS scenarios representa- 
tive of smaller and smaller incremental shifts in growth 
and size at maturity were added to determine the level 
of sensitivity (the breaking point), and those new growth 
matrices and sizes at maturity from these analyses 
were used in the UMM. For example, if historical fish- 
ery statuses from the UMM in which outputs from IBLS 


Table 2 


Settings and data that were consistent across scenarios used in the University of Maine Lobster Stock 
Assessment Model for American lobster (Homarus americanus) in the Gulf of Maine. IBLS=individual- 
based lobster simulator model; SSB/R=spawning stock biomass divided by the number of recruits to the 
stock; MEDMR=Maine Department of Marine Resources; MADMF=Massachusetts Division of Marine 
Fisheries; NHFGD=New Hampshire Fish and Game Department; NEFSC=NOAA Northeast Fisheries 


Science Center. 
Setting or data category 


Period for time series 
No. of seasons 

No. of sexes 

Size range 

Size bin length 

Initial conditions 
Recruitment size 
SSB/R relationship 

No. of commercial fleets 
Commercial fleet selectivity at size 
Survey data sources 


Survey selectivity at size 
Fishing mortality rate 
Natural mortality rate 


Description 


1984-2013 

4 (each a 3-month time block, same as IBLS) 

1 (data averaged across males and females) 

53-223 mm carapace length 

5 mm carapace length 

First-year size composition from survey data 

53-73 mm carapace length 

None 

1 

Double logistic 

MEDMR Ventless Trap Survey (2006-2012) 

Spring MEDMR/NHFGD inshore bottom-trawl survey (2001-2013) 
Fall MEDMR/NHFGD inshore bottom-traw] survey (2000-2013) 
Spring MADMF bottom-trawl survey (1984-2013) 

Fall MADMF bottomtrawl survey (1984-2013) 

Spring NEFSC bottomtrawl survey (1984—2013) 

Fall NEFSC bottom-trawl survey (1984-2013) 

Double logistic 

Instantaneous 

0.15 year | 
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scenario 6 were used did not differ from the status from 
the UMM base case but the statuses from the UMM in 
which IBLS scenario 7 outputs were used did differ from 
the status from the UMM base case, the breaking point of 
sensitivity would lie somewhere between the shifts repre- 
sented in the 2 IBLS scenarios. 

The next step was to produce growth matrices and esti- 
mate a size at maturity for a shift representative of half- 
way between these 2 shifts. For molting probability, the 
shift in this step was the average probability of both G1 
and G2 for each season since the last molt. For molt incre- 
ment probability, this shift was the average probability 
of both G1 and G2 for each size increase in millimeters. 
For size at maturity, it was simply the average of 85.21 
and 79.61 mm CL. Historical fishery statuses and BRPs 
were then calculated for these new UMM scenarios. Ret- 
rospective patterns were also evaluated, and results from 
these tests can be found in Supplementary Figures 2-12 
(online only). To further determine the breaking point of sen- 
sitivity, we used these new UMM scenarios in the methods 
described in the previous paragraph in place of either the 
UMM scenario in which IBLS scenario 6 was used or the 
UMM scenario in which IBLS scenario 7 was used (depend- 
ing on whether the results of these new scenarios were sig- 
nificantly different from the UMM base case). This process 
of determining a breaking point of sensitivity was repeated 
and continued until a breaking point within one-sixteenth 
of a shift was found. The results of the sensitivity analyses 
described in this section can help determine the sensitivity 
of the UMM to growth and size at maturity and can help to 
focus future research efforts toward direct linkages of cli- 
mate change to these life history parameters for use in stock 
assessment of American lobster. 


Results 


Target and threshold BRPs for the 7 UMM scenarios 
can be found in Table 3, and all accompanying reference 
abundance plots showcasing historical fishery statuses 
as compared to those of the UMM base case can be found 


in Figure 2. Terminal-year stock statuses did not change 
across any of the 7 UMM runs in this study (Table 4). Most 
alterations in historical reference abundance from the 
UMM base case appeared to be related to magnitude: con- 
sistent overestimations of abundance per year (except in 
the UMM scenario in which outputs from IBLS scenario 2 
were used, where abundance was underestimated) but 
similar temporal trends, with slight alterations causing 
some discrepancies in historical fishery statuses. There 
was no clear pattern in how each instance of discrep- 
ancy across the runs either improved or worsened stock 
status estimates. Instances of consecutive years differing 
from the UMM base case are much more relevant because 
these instances are indicative of larger trend-based dif- 
ferences and not simply 1-year lags that seem to be the 
reason behind solitary instances of differing years. These 
instances of differences across consecutive years appeared 
in only one UMM scenario: scenario 7. In this scenario, a 
growth shift of G2 and a size at maturity of 79.61 mm CL 
were used (the large climate effect scenario). 

Given that a change in size at maturity of over 10 mm CL 
did not appear to cause differences between consecutive 
years in reference abundance independent of a change in 
growth, a sensitivity analysis was not conducted for this 
variable. Likewise, changes in growth independent of 
size at maturity did not appear to cause consecutive-year 
differences. Therefore, a sensitivity analysis was not con- 
ducted for growth independent of size at maturity. 

In the biologically realistic scenarios (UMM scenarios in 
which data from IBLS scenarios 6 and 7 were used), the 
combination of G1 and the size at maturity of 85.21 mm 
CL resulted in no consecutive-year differences in histori- 
cal fishery status when compared to use of the UMM base 
case. However, the combination of G2 and the size at 
maturity of 79.61 mm CL resulted in consecutive-year dif- 
ferences compared to use of the UMM base case. There- 
fore, the breaking points of sensitivity existed somewhere 
between a small climate effect scenario (G1 and the size at 
maturity of 85.21 mm CL) and a large climate effect sce- 
nario (G2 and the size at maturity of 79.61 mm CL). 
Results from this sensitivity analysis for these biologically 


Table 3 


Target and threshold biological reference points (BRPs) for abundance (in millions 
of individuals) for all University of Maine Lobster Stock Assessment Model (UMM) 
scenarios for American lobster (Homarus americanus) in the Gulf of Maine. In each 
UMM scenario, the growth transition matrices and size at maturity produced with 
the corresponding scenario of the individual-based lobster simulator model were 
used (Table 1). The time series used in the UMM is from the period from 1984 
through 2013. 


UMM scenario 


BRP 1 2 3 4 5 6 7 


2857.5 
2170.5 


1403.6 1067.0 
1074.9 788.9 


Target 976.0 823.1 1256.3 1078.1 
Threshold 707.8 579.6 919.6 801.2 


246 


Scenario 2 


Reference abundance 
500 1000 


0 


Scenario 4 


® 
1S) 
= 
ia) 
xo) 
= 
= 
2 
ca) 
® 
oO 
ion 
® 
—_ 
® 
— 
) 
o 


Scenario 6 


Reference abundance 
500 1000 


0 


2013 


Reference abundance Reference abundance 


Reference abundance 


Fishery Bulletin 120(3-—4) 


Scenario 3 


Scenario 5 


Scenario 7 


S 
i=) 
2) 
= 
co) 

1 


984 


Figure 2 


Estimated reference abundance (in millions of individuals) of American lobster (Homarus americanus) in 
the Gulf of Maine during 1984—2013 for 6 scenarios in the University of Maine Lobster Stock Assessment 
Model (UMM) corresponding to scenarios in Tables 1 and 3 (black lines). Estimated abundance under 
the first scenario or the UMM base case is provided for comparison (dotted lines). The gray bars indicate 
years for which the historical fishery status (as calculated from the biological reference points in Table 3) 
is different from that of the UMM base case for the same year. Note the differences between panels in the 


range of values on the y-axis. 


realistic scenarios can be found in Figure 3. The final 
breaking point was between the growth shifts of 1.4375 
and 1.5000 and the size at maturity values of 82.41 and 
82.76 mm CL. 


Discussion 


Traditional sensitivity analyses are bottom-up: they are 
designed to determine how model output changes when 
specific parameters are altered (Booshehrian et al., 2012; 
Salciccioli et al., 2016). This practice is commonly used in 


stock assessment procedures to determine model stability 
and quantify uncertainty (Rosenberg and Restrepo, 1994; 
Hilborn, 2003; Salciccioli et al., 2016). The UMM scenar- 
ios in this study that involved IBLS scenarios 1—5 are an 
example of this classic type of analysis. The UMM scenar- 
ios in this study that involved IBLS scenarios 6—7, however, 
represent a top-down approach to sensitivity analysis. In 
this top-down approach, a large model-free mechanism con- 
trolled how multiple variables changed together and would 
affect model results. In this type of approach, the aim is to 
answer the question of how sensitive the model is to this 
large mechanism, which was climate change in this study. 
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Table 4 


Terminal-year stock abundance (in millions of individuals) and stock status from each University of 
Maine Lobster Stock Assessment Model (UMM) scenario for American lobster (Homarus americanus) 
in the Gulf of Maine. Also presented are proportions of each estimate of terminal-year abundance in 
relation to the target reference point. Stock status is presented in comparison to the biological reference 
points in Table 3. In each UMM scenario, the growth transition matrices and size at maturity produced 
with the corresponding individual-based lobster simulator model were used (Table 1). The time series 


used in the UMM is from the period from 1984 through 2013. 


Terminal-year 
parameter 1 Zi 


Abundance estimate 1594.5 1217.9 
Stock status >target >target 
Abundance proportion 1.63 1.48 


Climate change affects molting probability, molt incre- 
ment probability, and size at maturity of American lobster 
together. Therefore, this type of analysis is important to 
determine these cumulative effects on model outputs, suc- 
ceeding where traditional sensitivity analyses fail. This 
type of analysis is sometimes referred to as a global sen- 
sitivity analysis and is very rarely used in fishery stock 
assessments (Lehuta et al., 2010; Saltelli et al., 2019; 
Garcia, 2020). We agree with the notion of Saltelli et al. 
(2019) that a lack of the use of this method throughout 
the fields of environmental science and biology is concern- 
ing. We further postulate that both a bottom-up approach 
and a top-down approach may be beneficial and increas- 
ingly imperative in a changing world to ensure that a 
stock assessment model is stable under ensemble changes 
brought by large mechanisms. 

The sensitivity of the UMM to growth and size at matu- 
rity is relatively and biologically low. Size-at-maturity 
values associated with breaking points in the biologically 
realistic scenarios are not expected to reach such low lev- 
els for at least 50 years (Le Bris et al., 2017; IPCC, 2019; 
Brickman et al., 2021). The relationship of growth of 
American lobster to temperature and climate change are 
well-known (Aiken, 1977; Le Bris et al., 2017), but strict 
predictions cannot be easily extrapolated and may be less 
appropriate (Punt et al., 2014). This information, coupled 
with the fact that most information on these parameters 
found in laboratory settings may not be directly applicable 
to scenarios designed to represent conditions in the wild 
(Jury and Watson, 2013), means that forecasting growth 
and size at maturity of American lobster is incredibly 
challenging. An advantage of our modeling framework is 
that strict relationships of tested parameters to the large 
mechanism (e.g., climate change or temperature change) 
are not necessary. The framework is not meant to be 
used to determine future changes to modeling efforts, but 
rather results from using it can highlight the limitations 
of a stock assessment under climate change. 

It is important to consider that our study was focused on 
abundance as an output of interest because it is used to set 


1687.7 
>target 
1.34 1.64 1.61 1.50 1.33 


UMM scenario 


4 5 6 7 


3792.7 
>target 


1599.5 
>target 


1769.5 
>target 


2265.5 
>target 


management targets and thresholds. However, climate- 
driven changes to maturity and growth likely affect other 
aspects of a population and fishery dynamics, such as 
spawning stock biomass, recruitment, and fishing mortal- 
ity. The aim of our study was not to analyze the effects of cli- 
mate change on the full suite of outputs for population and 
fishery dynamics but rather was to understand whether 
determinations of fishery stock status would have been dif- 
ferent under the climate effects we simulated. 

Estimates of terminal-year stock status, most relevant 
to management of American lobster, did not change over 
all UMM scenarios in this study, indicating the robust- 
ness of the UMM to changes in life history parameters. 
However, results from using the combination of shifts in 
growth and size at maturity indicate differences in hind- 
casted fishery statuses. Consequently, scenarios tested in 
this study may not alter input data enough to produce dif- 
ferent results for current management, but given that his- 
torical deviations were present, caution should be given to 
the assumptions of low sensitivity. Deviations of historical 
stock statuses were mostly related to magnitude, repre- 
senting overestimations of abundance of American lobster 
throughout the time series but having very similar tempo- 
ral trends. This result is due to the use of relative BRPs 
calculated for each scenario as opposed to the use of static 
values over all scenarios. Lobster management, like much 
of fisheries management in general, is more concerned 
with trends (ASMFC, 2015) than absolute values. There- 
fore, large shifts in growth and size at maturity can alter 
model results but would not have severely affected histor- 
ical management decisions based on stock status. 

As expected, these UMM scenarios all had worse model 
fits (higher objective function values) than the UMM 
base case (see Supplementary Table) (online only). These 
worse fits are most likely due to the approximation of 
biologically unrealistic, freely estimated parameters 
in an attempt to fit to the input data while also using 
the growth and size-at-maturity data provided (Slezak 
et al., 2010). These differences in fit are not relatively 
high, even for the largest shifts in this study, but in work 
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Figure 3 


Estimated reference abundance (in millions of individuals) of American lobster (Homarus americanus) in the 
Gulf of Maine during 1984-2013 for the biologically realistic University of Maine Lobster Stock Assessment 
Model (UMM) scenarios used in the growth sensitivity analyses (black lines). The top-left panel represents a 
scenario between scenarios 6 and 7. Each of the other panels, the top-right panel, then the bottom panels from 
left to right, represents a consecutive scenario in the sensitivity analysis based on the scenario in the previous 


_ panel. In addition to the size at maturity (e.g., 82.41 cm carapace length [CL]), each panel title includes the shift 


from the individual-based lobster simulator model base case used to produce the abundance estimates shown, 
with the shift represented as a proportion between a full paired shift of molting probability by 1 year and molt 
increment probability by 1 size bin of 1 cm CL (G1) and a full paired shift of molting probability by 2 seasons and 
molt increment probability by 2 size bins (G2). Estimated abundance under the UMM base case (with original 
growth data [OG] and size at maturity) is provided for comparison in each panel (dotted lines). The gray bars 
indicate years for which the historical fishery status (as calculated from biological reference points in Table 3) is 
different from that of the UMM base case for the same year. The bottom panel shows the location of the breaking 
point of sensitivity within one-sixteenth of a shift (gray bar) in relation to the UMM base case data. In this bottom 
panel, the vertical lines represent the tests of partial shifts shown in the previous panels. Note the differences 


between panels in the range of values on the y-axis. 


on other models, this phenomenon should be considered. 
Caution should be used when using this approach, and 
careful attention should be paid to the freely estimated 
parameters of the model. 

It is important to note the combined effects of shifts in 
growth and size at maturity. The largest alteration in com- 
parison to the UMM base case occurred when the largest 
effects from growth and size at maturity were combined. 
In contrast, results from use of shifts that were smaller 
than that alteration indicate that combined effects may 
not be strictly additive, and future work should be focused 
on the complex relationships of growth, size at maturity, 
and temperature, especially as they pertain to the Ameri- 
can lobster stock assessment model. Quantifying the rela- 
tionships between these parameters and thermal habitat 
is a research priority (ASMFC, 2015), but another prior- 
ity is to develop modeling capacity to handle temporally 
dynamic life history parameters. If climate change affects 


key life history characteristics, traditional stock assess- 
ments in which static values for variables, such as growth, 
size at maturity, and others, are used may be inaccurate 
(Correa et al., 2021). 

A problem, however, highlights a key limitation of the 
analysis used in the current stock assessment model: cli- 
mate change affects many other aspects of the life history 
of American lobster besides growth and size at maturity. 
Key life history information, such as natural mortality 
(Mills et al., 2013) and recruitment (Goode et al., 2019; 
Tanaka et al., 2019), as well as important stock assess- 
ment model parameters, such as fishery-dependent and 
fishery-independent catchability or exploitation (Maunder 
et al., 2006; Conn, 2010; Shelton et al., 2014), are likely 
affected by changes in thermal habitat. Complex relation- 
ships with all of these factors exist, and a fully compre- 
hensive top-down sensitivity analysis would see climate 
effects on all of them together. 


Hodgdon et al.: Climate-driven changes in growth and size at maturity of Gulf of Maine lobster stocks 249 


The top-down approach used in the framework of our 
study can be used to examine a specific category of life 
history traits, but in future analyses, this work should 
be expanded to include effects from other life history and 
fishery components. Such a fully comprehensive anal- 
ysis, however, would be possible only after individual 
studies, such as the one we describe herein or was done 
by Hodgdon et al. (2020), in which the UMM was tested 
in the presence of thermal effects on survey catchability, 
are completed. Temporally dynamic life histories in stock 
assessment may require quantification of relationships 
with environmental conditions, but use of them would 
ultimately increase accuracy in model results and pre- 
cision of forecasts. Another avenue for future research 
would be the application of a management strategy eval- 
uation within the current framework. In this work, the 
IBLS would be used as an operating model so that results 
from the UMM could have a “true” population to use in 
comparisons. 

Ultimately, knowing the breaking points does not aid 
management if there is a lack of knowledge on the life 
history parameters tested a priori, specifically a lack 
of information concerning the relationship with each of 
them to thermal habitat and hypotheses as to the pre- 
dicted scale of future change. Foremost, the need to 
quantify the relationship that parameters related to life 
history traits of American lobster have with a changing 
climate is critical and is a concern that management 
shares (ASMFC, 2015). Such work is considered essen- 
tial because comparison of predicted changes to the stock 
assessment model’s breaking points aids in determining 
research necessity. If the breaking points are higher than 
the predicted changes, changes under climate change 
may not significantly affect stock assessments if data for 
the life history parameters (e.g., growth and size at matu- 
rity) are not updated. If the breaking points are lower 
than the predicted changes, modeling efforts with param- 
eter data that have not been updated may no longer yield 
accurate results, and the aim of future research should be 
to improve understanding of those parameters for which 
information is old. 


Conclusions 


In this study, we used a novel top-down sensitivity analy- 
sis for American lobster in the Gulf of Maine. This frame- 
work has great potential to better the management of 
the fishery for this species, but it should be expanded 
in future work. Lobster and, by extension, crustacean 
physiology and life history are directly linked to the envi- 
ronment and most often are consequences of thermal 
habitat (Madeira et al., 2012). As climate change alters 
the thermal habitats of crustaceans, the data used in 
stock assessment methods that rely on life history char- 
acteristics can become out of date. This issue can be mit- 
igated with persistent monitoring efforts and scientific 
research. However, many crustacean fisheries, even in 
well-funded jurisdictions such as the United States, have 


limited resources for such cost-intensive research efforts. 
The framework proposed in this paper can be used to 
possibly lessen research loads by prioritizing the input 
parameters to which a specific stock assessment model is 
most sensitive under the top-down mechanism of climate 
change. A complete analysis of dependent and indepen- 
dent effects from all variables together under this frame- 
work has the potential to aid management practices, to 
advance crustacean stock assessment, and to steer future 
research projects. 


Resumen 


Las especies de crustaceos son socioeconémica y ecolégi- 
camente cruciales en todo el mundo. Para los crusta- 
ceos, como ectotermos, el cambio climatico antropogénico 
amenaza con alterar significativamente las caracteristicas 
clave de su historia de vida, como la talla de madurez y 
la tasa de crecimiento. Dado que es dificil determinar la 
edad de los crustaceos, los datos de longitud se utilizan 
en las evaluaciones de sus poblaciones; sin embargo, los 
cambios en la maduracion y el crecimiento, por efecto del 
clima pueden influir en gran medida en el desempeno de 
los modelos estructurados por talla en la evaluacién de 
las poblaciones. Acoplamos modelos individuales y estruc- 
turados por talla para la langosta americana (Homarus 
americanus) frente al noreste de Norteamérica, en el golfo 
de Maine, para realizar un novedoso andalisis de sensib- 
ilidad del efecto de los parametros de entrada de madu- 
rez y crecimiento en los resultados del modelo. Para este 
analisis, utilizamos un enfoque ascendente (con paramet- 
ros modificados de forma independiente) y un enfoque 
descendente (con parametros modificados de forma con- 
junta a medida que se preveia su influencia por el cambio 
climatico). Encontramos que nuestro modelo de evalu- 
acién de la poblacién de langosta americana es resiliente 
a los cambios relativamente extremos de los parametros 
biol6gicos de entrada. Para la evaluacion de poblaciones de 
crustaceos con modelos estructurados por talla, recomen- 
damos ampliar los andlisis de sensibilidad para incluir la 
evaluacion de la influencia de los cambios provocados por 
el clima sobre los parametros de entrada basados en los 
rasgos del ciclo de vida que varian con el tiempo. 
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Abstract—Abundance indices from 
fishery-independent surveys are pre- 
ferred in stock assessments for their 
robust scientific designs that minimize 
uncertainty and bias. When sampling 
does not adhere to the design, research- 
ers employ techniques such as impu- 
tation or standardization to improve 
accuracy and reduce bias. We examined 
2 methods for adjusting for incomplete 
sampling within the Coastal Trawl 
Survey (CTS) of the Southeast Area 
Monitoring and Assessment Program— 
South Atlantic for 3 species commonly 
encountered in survey sampling, the 
Atlantic croaker (Micropogonias undula- 
tus), bluefish (Pomatomus saltatrix), and 
white shrimp (Litopenaeus setiferus): 
design-based imputation of missing 
data and standardization through the 
delta-generalized-linear-model approach. 
Additionally, we determined the effect 
of modifying the seasonal component of 
the survey design through retrospective 
simulation. For all 3 species, standard- 
ization improved precision in annual 
abundance estimates relative to val- 
ues estimated with the design-based 
method. When a stratum missed in sam- 
pling overlapped with an area or time 
of high variability for a species (e.g., 
2019), standardization did not improve 
precision over the design-based method. 
Results from examination of the effects 
of dropping entire seasons, because of 
funding or logistical challenges, indicate 
that rotating which season is dropped 
was the best approach to balancing 
characteristics of each species. Over- 
all, we recommend the standardization 
approach for accounting for missing 
data within the CTS time series. 
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Generally, the purpose of fishery- 
independent surveys is to provide pop- 
ulation data to stock assessments to 
reduce uncertainty in stock statuses 
and management measures informed 
by the stock assessments (Walters and 
Pearse, 1996). Reduction in uncertainty 
relies on fishery-independent surveys 
having robust scientific designs that 
are not influenced by management 
or fishing practices (Cochran, 1977; 
Williams and Carmichael’). To reduce 
assessment uncertainty, considerable 
effort has been applied to developing 
best practices for minimizing error in 
abundance indices derived from sur- 
vey data (Walters, 2003; Maunder and 
Punt, 2004; Shelton et al., 2014). Addi- 
tionally, much work has been done 
to address bias in index calculation 
due to spatial variability in sampling 
or fishing and to ensure that effort 
is defined appropriately (Campbell, 
2004, 2015). Precision can be gained in 


1 Williams, E. H., and J. Carmichael (eds.). 
2009. South Atlantic fishery indepen- 
dent monitoring program workshop final 
report, 85 p. South Atl. Fish. Manag. 
Counc. and Natl. Mar. Fish. Serv., South- 
east Fish. Sci. Cent., Beaufort, NC. 
[Available from Southeast Fish. Sci. 
Cent., Natl. Mar. Fish. Serv., 101 Pivers 
Island Rd., Beaufort, NC 28516.] 


fishery-independent surveys either by 
increasing sampling effort or by using 
a stratified-random sampling design 
to optimize effort allocation (Xu et al., 
2015). However, the realities of fund- 
ing, weather, and vessel reliability 
and availability more often than not 
result in decreased or incomplete 
sampling effort decreasing precision 
of a survey. In particular, changing 
environmental conditions and fund- 
ing concerns that affect completion of 
surveys are highly pervasive issues 
among surveys, and many survey 
programs are facing hard decisions 
regarding reducing effort while still 
providing comparable time series. In 
the absence of increased precision 
from a priori design, researchers may 
turn to analytical tools. 

Stock assessment models often rely 
on the inclusion of catch rate time 
series (either fishery-independent or 
fishery-dependent) that are propor- 
tional to population abundance (gener- 
ally referred to as indices of abundance) 
(Francis, 2011). Often it is assumed 
that these time series are proportional 
to stock size (Hilborn and Walters, 
1992; Arreguin-Saénchez, 1996; Quinn 
and Deriso, 1999), on the basis of the 
sampling design, stability of catchabil- 
ity, and consistency in trends among 
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indices of abundance (Quinn and Deriso, 1999). For an 
index of abundance to reflect stock size trends, the data 
must be collected according to a robust scientific design 
consistently through time. When deviations from the 
design occur, those deviations should be accounted for in 
some way. Standardizing catch data increases the accu- 
racy of estimates (Walters, 2003; Maunder and Punt, 2004) 
and allows accounting for incomplete sampling during the 
time series. Typical statistical standardization techniques 
(e.g., generalized linear models or generalized additive 
models) provide the opportunity to quantify the effects of 


South Carolina 


81°W 79°W 


either sampling or environmental characteristics on esti- 
mates of abundance and to essentially correct those esti- 
mates for deviations in sampling or environmental 
conditions (Wilberg et al., 2009). 

Off the Atlantic coast of the southeastern United 
States, fish, elasmobranch, and invertebrate species in 
nearshore, coastal waters are surveyed through fishery- 
independent sampling of the Southeast Area Monitor- 
ing and Assessment Program—South Atlantic’s Coastal 
Trawl Survey (SEAMAP-SA CTS). Sampling of trawlable 
habitats through the CTS began in this region (Fig. 1) in 
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Figure 1 


Map of study area showing the stratum boundaries and locations of stations of the 
Coastal Trawl Survey of the Southeast Area Monitoring and Assessment Program— 
South Atlantic where Atlantic croaker (Micropogonias undulatus), bluefish (Pomatomus 
saltatrix), and white shrimp (Litopenaeus setiferus) were sampled from 1990 through 
2019. Numbers are those assigned to strata. Stations are located in trawlable habitat 
from 4.6 to 9.1 m off the Atlantic coast of the southeastern United States. 
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1986, and since 1990 the CTS has monitored abundance 
trends in shallow depths (<9.1 m) using a standardized 
stratified-random sampling design. Each year, the num- 
ber of stations allocated for sampling is determined by 
available funding, but the spatial and temporal portions 
of the survey design have remained consistent over time. 
This consistency makes it possible for estimates based on 
catch rate time series from this survey to serve as indices 
of abundance for stock assessments. 

Historically, the rates at which stations were missed 
were low, and the missed stations were generally spread 
out among strata instead of entire strata being missed. 
More recently, however, the CTS has experienced chal- 
lenges to sampling all of the stations allocated for sam- 
pling as a result of the increased frequency of weather 
conditions unsuitable for sampling (including prolonged 
periods of winds and sea states above sampling cut-offs 
and major hurricanes), damaged gear (caused by debris 
from river run-off following hurricanes and major rain 
events), logistical challenges (including low tidal ampli- 
tude preventing departures), and mechanical failures 
making the aging vessel unavailable for sampling. In 
some seasons during 2018 and 2019, these challenges 
meant that entire strata were not sampled. Given these 
challenges, it has become more likely that annual abun- 
dance estimates based on survey data from years with 
completely unsampled strata or high numbers of unsam- 
pled stations will not be proportional to true abundance 
or will not index the same portion of the population as 
estimates based on data from other years with more com- 
plete sampling of strata. 

Herein, we present catch and effort time series for 
3 species of management interest commonly encountered 
in CTS sampling: the Atlantic croaker (Micropogonias 
undulatus), bluefish (Pomatomus saltatrix), and white 
shrimp (Litopenaeus setiferus). These species have a 
variety of distributions (spatially and temporally) (senior 
author, unpubl. data), and incomplete sampling may 
affect abundance estimates in different ways. In light of 
recent difficulties completing CTS sampling, we present 
annual abundance for these species in 2 ways, a nominal, 
design-based estimate and a standardized estimate, to 
assess if these time series need correcting through stan- 
dardization because of incomplete sampling. In addition, 
we provide results of our examination of the effects of 
sampling and environmental covariates on the standard- 
ized estimate of abundance to identify potential drivers 
of abundance and the need for correction for each spe- 
cies. Because the challenges for completing survey sam- 
pling as designed are expected to continue, changes to 
the survey design, such as dropping a sampling season, 
are under consideration. To determine the effects of these 
potential design modifications on abundance estimates, 
we also conducted retrospective simulations on the stan- 
dardized indices of abundance for all 3 species. On the 
basis of these analyses, we provide a recommendation on 
which estimate is most appropriate for use as an index 
of abundance and on whether there is a preferred design 
modification strategy. 


Materials and methods 
Study design and gear 


Sampling for the CTS, a fishery-independent research 
program, was conducted by the South Carolina Depart- 
ment of Natural Resources (SCDNR) in coastal Atlantic 
waters off the southeastern United States between Cape 
Hatteras, North Carolina, and Cape Canaveral, Florida 
(Fig. 1). Consistent sampling methods were in place from 
1990 through 2019 and were used to target the full spa- 
tial range in each of 3 seasons: spring (April and May), 
summer (July and August), and fall (September and 
November). Trawl surveys were conducted at randomly 
selected stations from a pool of stations within 24 strata 
based on latitude and depth (Fig. 1) between the inshore 
4.6-m depth contour and the offshore 9.1-m depth contour. 
To reduce variability of the data, the method of allocat- 
ing stations was changed in 2001 from proportional allo- 
cation (based on the total surface area of each stratum) 
to optimal allocation (Thompson, 1992), with higher effort 
allocated to strata with historically higher variability and 
with the number of stations allocated within each stra- 
tum determined anew each year. The total number of allo- 
cated stations per season has ranged between 78 and 112, 
depending on funding and other survey priorities, but the 
spatial footprint of the survey has remained consistent. 

Sampling was conducted during daylight hours 
(between 1 h after sunrise to 1 h before sunset) on board 
the R/V Lady Lisa, a 22.9-m wooden-hulled, double- 
rigged St. Augustine shrimp trawler owned and oper- 
ated by the SCDNR Marine Resources Division, by using 
a pair of 22.9-m mongoose-type Falcon trawl! nets (Beau- 
fort Marine Supply Inc.”, Beaufort, SC) without turtle 
excluder devices (Willis et al., 2015; Zimney’). The body 
of the trawl net was constructed of #15 net twine with 
47-mm stretch mesh, and the codend was constructed of 
#30 net twine with 41-mm stretch mesh. At each station, — 
the pair of nets were towed for 20 min, excluding wire- 
out and haul-back times, with a target speed of 1.3 m/s, 
relative to the bottom. 

The catch from each net was processed independently 
and assigned a unique collection number. The contents of 
each net were sorted to species (with limited exceptions 
sorted to genus or family only), and the total biomass and 
number of individuals were recorded. When trawl nets 
contained high volume catches, selected species were 
removed (e.g., endangered species and species that posed 
a risk to staff during handling), the remaining net con- 
tents were placed into shrimp baskets and weighed, and 
a randomly selected basket was sorted and processed as 
described above. Abundance and biomass data for each 


2 Mention of trade names or commercial companies is for identi- 
fication purposes only and does not imply endorsement by the 
National Marine Fisheries Service, NOAA. 

3 Zimney, A. 2021. SEAMAP-SA Coastal Trawl Survey data and 
sample collection methods. Southeast Data, Assessment, and 
Review SEDAR78-WP01, 4 p. [Available from website.] 
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species were then used to estimate the total abundance 
and biomass of each species in the full catch by using the 
ratio of the weight of subsampled catch placed in baskets 
to the total weight of the full catch. When large numbers of 
an individual species occurred in the processed catch, all 
individuals of that species were weighed, and then a hap- 
hazardly selected subsample was counted (subsamples 
were approximately 30-60 individuals). The total number 
of individuals of a species in the catch was then estimated 
by using the ratio of the weight of the processed species 
subsample to the total species weight. 

Hydrographic data, measurements of surface and bot- 
tom temperature and salinity, were logged at each sta- 
tion, with a Van Dorn water sampler (EKijkelkamp North 
America Inc., Wilmington, NC) in 1990-1992, an SBE 19 
SeaCAT Profiler CTD in 1993-2005 and 2017 (Sea-Bird 
Scientific, Bellevue, WA), or an SBE 19plus SeaCAT 
Profiler CTD in 2006—2019 (V1 or V2, Sea-Bird Scientific). 
From these CTD casts, we extracted the deepest measure- 
ment of temperature and salinity, as long as it was within 
1 m of the bottom, to use in analyses (hereafter, these mea- 
surements are referred to as bottom temperature and bot- 
tom salinity). 


Data and nominal abundance estimation 


The data available for use in abundance estimation 
included a unique collection number, date of deployment, 
season, tow duration (in minutes), stratum, station, lati- 
tude, longitude, depth, bottom temperature, bottom salin- 
ity, and the number of individuals and aggregate weight 
of each species captured. Data from the paired nets were 
pooled for analysis to form a standard unit of effort at 
each station sampled (a trawl tow). Estimates of abun- 
dance were expressed as the number of individuals per 
hectare. Estimated area swept by a net was calculated by 
multiplying the expected average width of the net open- 
ing (13.5 m), based on Stender and Barans (1994), by the 
distance in meters trawled and dividing the product by 
10,000 m*/ha. Then, the estimated area swept for each 
of the paired nets was combined. If area swept could not 
be accurately estimated (e.g., if a U-turn was executed to 
avoid entanglement with a gill net), data for that tow were 
omitted from analyses. 

From 1990 through 2019, sampling was conducted at 
8403 of the 8568 allocated stations, and 24 strata were 
fully missed (Table 1, Fig. 1): 1 stratum (21) in fall 1990, 
2 strata (65 and 67) in spring 2013, 4 strata (37, 39, 65, and 
67) in spring 2018, 7 strata (29, 57, 59, 61, 68, 65, and 67) 
in fall 2018, 3 strata (63, 65, and 67) in spring 2019, and 
7 strata (21, 23, 25, 27, 63, 65, and 67) in fall 2019. Where 
entire strata were missed, the average long-term nominal 
abundance was calculated for each species in each stratum 
and season and imputed for the missing value as the nom- 
inal abundance. In total, 58 trawl tows were excluded from 
analyses (Table 1). Nine tows were excluded because the 
area swept could not be accurately estimated, and 14 tows 
were eliminated because the tow duration was not equal 
to 20 min, indicating deviation from normal operations. 


Additionally, 35 tows were excluded from standardization 
because covariate information was missing. 

Annual nominal abundance (A) was calculated by deter- 
mining the sum of the number of individuals (individuals) 
caught per hectare (area swept) and dividing it by the total 
number of tows (¢) for a given combination of stratum (st), 
season (se), and year (y): 


A, =Yia 


individuals, «tse y 
= ie . (1) 
? J 
area SWEPLt, «t, se,y 


The nominal index of abundance was then normalized by 
dividing the annual nominal abundance by the overall mean 
abundance for the full time series, producing a value of rel- 
ative abundance for each year. These values provide refer- 
ence points for individual years in relation to the time series, 
with a value of 1 being the long-term mean. 

In recent years (2015-2019), sampling for the CTS has 
not been completed at all of the allocated stations for the 
reasons described in the “Introduction.” Consequently, 
entire strata have been missed in some seasons primarily 
in the outer ranges of the study area where longer windows 
of suitable weather conditions are needed for travel and 
sampling (Table 1). Where an entire stratum was missed 
during a season, the long-term average nominal abundance 
was calculated for each stratum and season and imputed 
as a proxy for the nominal abundance of the missed stra- 
tum for use in Equation 1 (Little and Rubin, 1987; Walters, 
2003). The long-term average was used rather than the 
average for the most recent years because using only data 
from recent years could have introduced bias if a stratum 
was sampled unusually early or late in a season (timing 
can vary by up to 6 weeks) or was missed in several of the 
recent years, as occurred for strata 65 and 67, for example. 


Standardization with delta-generalized linear models 


Because of incomplete seasonal and regional sampling cov- 
erage, annual abundance was standardized among years 
through the “delta-GLM” method (Lo et al., 1992; Dick, 
2004; Ballenger et al.*), in which the standardized abun- 
dance is produced by using the product of predicted val- 
ues from 2 generalized linear models (GLMs). In the first 
GLM, species presence and absence is the response vari- 
able, sampling or environmental covariates are included, 
and the binomial error distribution is used. In the second 
GLM, only data from the tows for which the species of 
interest were present are used, and the number per hect- 
are is the response variable. Sampling and environmental 
covariates are also included in this second model, and this 
positive GLM is fitted with either the gamma or lognormal 
error distribution. Both gamma and lognormal error dis- 
tributions were considered for the positive GLM because 
preliminary analyses found that both error distributions 


* Ballenger, J.,T. Smart, K. Kolmos, and M. Reichert. 2013. Trends 
in relative abundance of gray triggerfish in waters off the SE 
US based on fishery-independent surveys. Southeast Data, 
Assessment, and Review SEDAR32-DW04, 77 p. [Available from 
website. ] 
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Table 1 


The number of stations allocated for sampling each year, the number of stations at which survey tows were completed, 
the survey strata with missed stations, and the number of stations with available data that were included in analy- 
sis for standardization of indices of abundance for Atlantic croaker (Micropogonias undulatus), bluefish (Pomatomus 
saltatrix), and white shrimp (Litopenaeus setiferus) sampled during the Coastal Trawl Survey of the Southeast Area 
Monitoring and Assessment Program—South Atlantic off the Atlantic coast of the southeastern United States from 
1990 through 2019. Percent positive rates, the percentage of tows with species present, are given for each species in 


each year. An asterisk (*) indicates that a stratum was missed in at least one season in the given year. 


No. of No. of Percent positive rate (%) 
No. of stations with stations 
stations complete with data Atlantic White 
Year allocated sampling Strata with missed stations used croaker Bluefish shrimp 
1990 234 2a) 21*, 61 227 65.2 43.2 50.2 
E99L 234 233 61 231 64.5 36.8 47.6 
1992 234 234 — 232 56.0 40.1 aie 
1993 234 234 - 233 50.2 30.0 39.9 
1994 234 234 ~ 233 53.2 32.6 43.8 
1995 234 234 _ 232, 59.1 47.8 51.7 
1996 234 234 - 229 58.5 45.3 48.3 
1997 234 234 232 43.8 30.5 47.2 
1998 234 234 = 232 64.4 33.5 57.9 
1999 234 234 - 233 44.6 35.6 55.8 
2000 234 234 _ 234 44.9 37.2 54.7 
2001 306 306 — 295 63.0 46.9 44.6 
2002 306 306 — 301 45.9 23.9 49.2 
2003 306 306 ~ 304 63.1 46.1 43.5 
2004 306 306 -- 304 54.9 36.9 46.7 
2005 306 306 - 301 58.5 28.8 39.5 
2006 306 306 _ 303 54.8 20.7 45.6 
2007 306 306 o 305 54.2 De 55.6 
2008 306 306 - 304 53.3 36.3 62.1 
2009 336 336 -- 335 60.9 40.9 51.6 
2010 336 336 _ 335 50.1 40.0 40.3 
2011 336 336 a 336 59.8 38.1 44.9 
2012 336 336 ~ 336 69.9 33.3 64.0 
2013 306 295 63, 65*, 67* 294 83.1 39.7 62.4 
2014 306 306 ~ 306 70.9 38.9 43.8 
2015 336 329 35, 63, 65 329 71.4 35.0 50.2 
2016 336 331 29, 33, 35, 39 330 70.0 31.5 64.2 
2017 306 294 21, 31, 35, 59, 61, 63, 67 293 71.8 22.1 65.6 
2018 306 228 21923, 25,7720", 31, 53,55, 225 78.9 Za.2 60.5 
3171397 4151, 53, 00,1", 
pot,O47 03°, 65°, 67" 
2019 306 258 7h Nga S hanes git bagta 29 es 2 lag os 258 84.9 24.0 69.0 
63*, 65*, 67* 
Total 8568 8403 8345 61.1 34.8 51.3 


fit data from the CTS for several species relatively well. 
Additionally, although other zero-inflated models have 
been used for standardization, the design of the CTS works 
well with the delta-GLM framework, and other models do 
not provide large improvements in variability for common 
species such as those examined here (Smart and Zimney’). 


° Smart, T,, and A. Zimney. 2021. Spanish mackerel indices of abun- 
dance in U.S. South Atlantic waters based on the SEAMAP-SA 
fishery-independent Coastal Trawl Survey. Southeast Data, 
Assessment, and Review SEDAR78-WP02, 22 p. [Available from 
website. ] 


Both error distributions were compared with the Akaike 
information criterion (AIC) in identical base models, and 
the one with the lowest AIC value was chosen as the model 
to use for analysis (Akaike, 1973). 

Year, season, stratum, bottom temperature, and bottom 
salinity were included as model covariates. Season was 
defined as spring, summer, and fall as described earlier. 
Stratum was previously defined in the first paragraph of 
the “Materials and methods” section and in Figure 1 and 
served as a proxy for location (latitude and longitude). 
Observed values of bottom temperature and bottom salin- 
ity were binned by using quantiles, with 5 bins for bottom 
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temperature (<20.2°C, 20.2—23.0°C, 23.1-26.5°C, and 
>26.6°C) and 2 bins for bottom salinity (<34.9 and 234.9). 
The data for these covariates were binned rather than left 
continuous for consistency with accepted index formula- 
tions for the CTS in recent stock assessments (Smart and 
Boylan®). The addition of these covariates also is consis- 
tent with recent standardization techniques used in many 
recent assessments in the region. The delta-GLM for all 
species started with the same base model: 


A, = Pry, se, st, temp, sal) x A,(y, se, st, temp, sal), (2) 


where Pr = the likelihood of a species being present in year y, 
season se, stratum st, bottom temperature 
temp, and bottom salinity sal; and 

A, = the abundance of a species when it is present 
in year y, season se, stratum st, bottom tem- 
perature temp, and bottom salinity sal. 


No interaction terms were included in the delta-GLM for- 
mula for consistency with recent accepted index formula- 
tions (Smart and Boylan®; Smart and Zimney”). Results of 
preliminary data analysis indicate no obvious interactive 
trends for any individual species examined in this study, 
and no multicollinearity was found for covariates included 
in the equation. 

The covariates included in the final delta-GLM models 
(both the binomial GLM and positive GLM) were chosen 
by using AIC with backward selection, with the excep- 
tion that year was always included in each model to allow 
an annual abundance value to be produced. The effect 
of covariates on standardized abundance estimates for 
each species was determined by using separate analy- 
ses of variance (ANOVA) for both the presence of a given 
species and the abundance of that same species esti- 
mated with model covariates as fixed factors. The final 
delta-GLM standardized abundance index for each spe- 
cies is the product of the mean responses from the linear 
predictors of year and any selected covariates from the 
2 GLMs. Jackknifing was used to determine coefficients 
of variation, standard error, and standard deviations for 
each delta-GLM analysis. Residual and quantile-quantile 
(Q-Q) plots were used as diagnostic tools to determine 
the fit of the models. All analyses were performed in R, 
vers. 4.1.1 (R Core Team, 2021). 

The delta-GLM used in this study was based primarily 
on code adapted from Dick (2004). As with the nominal 
index, the delta-GLM standardized index was normalized 
by dividing the annual standardized abundance by the 
overall mean standardized abundance for the time series, 
with a value of 1 being the mean for the time series. Nom- 
inal abundance calculations and delta-GLM analysis were 
performed for 3 species of management interest: the Atlan- 
tic croaker, bluefish, and white shrimp. Because these spe- 
cies also have generally been commonly caught in tows 


6 Smart, T. I., and J. Boylan. 2013. King mackerel index of abun- 
dance in coastal US South Atlantic waters based on a fishery- 
independent trawl survey. Southeast Data, Assessment, and 
Review SEDAR38-DW-11, 39 p. [Available from website.] 


over the full length (1990-2019) of the CTS time series, 
their distributions and abundances lend themselves to the 
analyses outlined herein. Over the period of the CTS time 
series, the Atlantic croaker was often the most abundant 
species caught during CTS sampling annually, having 
widespread spatial and temporal distributions. Bluefish 
primarily were caught in northern strata in spring and 
fall. Given the recent sampling challenges during these 
seasons, estimates of the abundance of bluefish have the 
potential to be heavily affected by incomplete sampling. 
The white shrimp was another widely distributed species, 
largely occurring from northern Florida to central North 
Carolina. 

In addition to comparison of design-based and stan- 
dardized indices, the effects of several potential survey 
design modifications were investigated for each species. 
Design modifications under consideration are primarily 
focused on dropping a season from sampling in any given 
year. Potential scenarios for this design change are per- 
manently dropping a season (3 scenarios: spring, summer, 
or fall), rotating the season that is dropped each year, or 
randomly selecting a season to drop each year. We retro- 
spectively applied these 5 scenarios to the raw data for 
the 3 species during the period 2014-2019 to investigate 
the effect of the design for years when sampling coverage 
was severely incomplete (i.e., in 2018 and 2019) or was 
complete or mostly complete (i.e., in 2014-2017). Once 
these scenarios were applied to the raw data, we ran the 
delta-GLM analysis again for the full time series and 
examined the standardized abundance estimates and 
variability for all 5 scenarios and for the current survey 
design that includes all seasons. 


Results 


All covariates (year, season, stratum, bottom temperature, 
and bottom salinity) were selected for inclusion by using 
the AIC for both the binomial and positive delta-GLMs for 
the Atlantic croaker, bluefish, and white shrimp. For all 
3 species the lognormal error structure had the lowest AIC 
value, indicating that it was the most appropriate error 
distribution to be used for the positive model. 


Atlantic croaker 


The Atlantic croaker was the most abundant species 
caught in CTS sampling, of the 3 species examined in this 
study. It was also the most common, occurring in 61% of 
all tows between 1990 and 2019 (Table 1). Atlantic 
croaker were distributed relatively evenly, both spatially 
and temporally, throughout the study area with speci- 
mens caught in all seasons and strata (Fig. 2, A and B). 
Atlantic croaker generally were more abundant in sum- 
mer than in spring and least abundant in fall. Although 
this species occurred in all strata, abundance of Atlantic 
croaker increased from the southern portion of the study 
area to the northern portion. Atlantic croaker were abun- 
dant in a range of bottom temperatures but were more 
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Figure 2 


Nominal abundance of Atlantic croaker (Micropogonias undulatus), by (A) season, (B) stratum, (C) bottom tempera- 
ture, and (D) bottom salinity, based on data from sampling off the Atlantic coast of the southeastern United States 
from 1990 through 2019, with design-based imputation of missing values. (E) Nominal and standardized relative 
abundances of Atlantic croaker estimated by using the design-based method and delta-generalized linear models, 
respectively, are provided by year. Error bars indicate the standard errors of the means. The horizontal gray line indi- 


cates the long-term mean of 1 for relative abundance over the full period of the time series. 


abundant in waters with bottom salinities greater than 
34.9 than in waters with other salinities (Fig. 2, C and D). 

From 1990 through 2010, nominal abundance of Atlantic 
croaker remained relatively low compared with the mean 
nominal abundance for the time series, with the normal- 
ized nominal abundance below the time series mean for 
most years (Fig. 2E). After 2010, nominal abundance was 
more variable and had a general trend of values increas- 
ing to levels above the time series mean (Fig. 2E). 

The trend in delta-GLM standardized abundance of 
Atlantic croaker was similar to that in nominal abun- 
dance, although the trend in standardized estimates 
was less variable (Fig. 2E, Suppl. Table [online only]). The 
exception to this similarity is that, for 2019, the pre- 
dicted abundance of Atlantic croaker in the standard- 
ized index was almost twice the abundance predicted 
in the nominal index. This large difference between the 
nominal and standardized abundance values in 2019 
coincides with the only occasion when the variability 


in the standardized estimate was higher than that of 
the nominal estimate for this species, indicating that 
variability drove the disagreement between the esti- 
mates produced with the 2 methods. In 2018 and 2019, 
limited sampling occurred in the northernmost strata 
where Atlantic croaker were generally more abundant 
than in other strata. The raw average nominal abun- 
dance in 2019 (excluding the imputed values for missed 
strata) was the highest in the history of the CTS (359.8 
individuals/ha) (senior author, unpubl. data), possi- 
bly contributing to the standardized estimate being 
higher than the nominal estimate. In contrast, the data 
imputed for missed strata for nominal abundance were 
from years with lower abundance. All covariates were 
significant in the binomial and positive delta-GLMs 
for predicting the presence and abundance of Atlantic 
croaker (Table 2). Results from the use of diagnostics, 
such as residual and Q-Q plots, indicate reasonable fits 
for both models (Suppl. Fig. 1) (online only). 
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Bluefish 


Bluefish were collected in 35% of all CTS tows from 1990 
through 2019 (Table 1). Temporally, bluefish were more abun- 
dant in spring and fall (Fig. 3A). The spatial distribution of 


Results from analysis of variance for presence and abundance of Atlantic croaker (Micropogonias undulatus), 
bluefish (Pomatomus saltatrix), and white shrimp (Litopenaeus setiferus) from standardization with binomial 


Table 2 


and positive delta-generalized linear models. Data used in models are from sampling of the Coastal Trawl 


Survey of the Southeast Area Monitoring and Assessment Program—South Atlantic conducted during 1990- 


2019 off the Atlantic coast of the southeastern United States. The significance level is 0.05. 


Covariate 


Atlantic croaker 
Binomial model 
Null 
Year 
Season 
Stratum 
Temperature 
Salinity 
Positive model 
Null 
Year 
Season 
Stratum 
Temperature 
Salinity 
Bluefish 
Binomial model 
Null 
Year 
Season 
Stratum 
Temperature 
Salinity 
Positive model 
Null 
Year 
Season 
Stratum 
Temperature 
Salinity 
White shrimp 
Binomial model 
Null 
Year 
Season 
Stratum 
Temperature 
Salinity 
Positive model 
Null 
Year 
Season 
Stratum 
Temperature 
Salinity 


df 


29 


23 


29 


23 


29 


23 


29 


23 


29 


23 


Deviance 


413.00 
498.71 
1188.01 
11.06 
63.89 


2590.40 
1293.70 
4554.10 
140.80 
103.00 


202.78 
242.50 
547.31 
166.85 

25.37 


330.53 
157.77 
1093.65 
197,89 
4.43 


251.60 
385.29 
1217.70 
90.63 
65.95 


1244.47 
1436.22 
1545.00 
46.54 
81.40 


Residual Residual 
df deviance F oa 

8344 11,159.10 
8315 10,746.10 <0.01 
8313 10,247.40 <0.01 
8290 9059.40 <0.01 
8287 9048.30 0.01 
8286 8984.40 <0.01 
5092 31,609.00 
5063 29,019.00 19.61 <0.01 
5061 27,725.00 142.02 <0.01 
5038 23,171.00 43.47 <0.01 
5035 23,030.00 10.31 <0.01 
5034 22,927.00 22.62 <0.01 
8344 10,788.80 
8315 10,586.00 <0.01 
8313 10,343.50 <0.01 
8290 9796.20 <0.01 
8287 9629.30 <0.01 
8286 9604.00 <0.01 
2906 6282.10 
2877 5951.50 Tae <0.01 
2875 5793.70 49.95 <0.01 
2852 4700.10 30.11 <0.01 
2849 4502.20 ave? <0.01 
2848 4497.80 2.80 0.09 
8344 11,562.80 
8315 11,311.20 <0.01 
8313 10,925.90 <0.01 
8290 9708.20 <0.01 
8287 9617.60 <0.01 
8286 9551.60 <0.01 
4282 17,992.00 
4253 16,748.00 13.29 <0.01 
4251 15,312.00 222.40 <0.01 
4228 13,767.00 20.80 <0.01 
4225 13,720.00 4.81 <0.01 
4224 13,639.00 25:21 <0.01 


bluefish was primarily in the northern portion of the study 
area (strata >55), with the highest abundances occurring 
off North Carolina (Fig. 3B). Bluefish were more abundant 
in waters with bottom temperatures below 20.1°C and in 
waters with bottom salinities less than 34.9 (Fig. 3, C and D). 
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Figure 3 


Nominal abundance of bluefish (Pomatomus saltatrix), by (A) season, (B) stratum, (C) bottom temperature, and 
(D) bottom salinity, based on data from sampling off the Atlantic coast of the southeastern United States from 1990 
through 2019, with design-based imputation of missing values. (E) Nominal and standardized relative abundances of 
bluefish estimated by using the design-based method and delta-generalized linear models, respectively, are provided 
by year. Error bars indicate the standard errors of the means. The horizontal gray line indicates the long-term mean 
of 1 for relative abundance over the full period of the time series. 


Nominal abundance of bluefish was generally stable 
near the time series mean without large positive or neg- 
ative swings. Additionally, variability in nominal abun- 
dance of bluefish tended to be low when abundance was 
below the mean and high when abundance was above 
the mean (2004—2013) (Fig. 3E). Nominal abundance 
was similar to or below the time series mean from 2014 
through 2019. 

The trend in delta-GLM standardized abundance 
of bluefish was similar to that in nominal abundance, 
although with less variability, particularly when com- 
pared with trends for the years of peak nominal abun- 
dance: 1995, 2004, 2005, 2009, 2010, and 20138 (Fig. 3E, 
Suppl. Table [online only]). One exception is the estimate 
for 2001, which is almost twice as high with higher vari- 
ability than that of the nominal estimate for that year. 
Although the error bars overlap, standardization also 
reduced the abundance estimates in 2009 and 2010 
relative to the nominal estimates for those years. All 


covariates were significant in both the binomial and pos- 
itive GLMs, except for bottom salinity, which was not sig- 
nificant in the positive model (Table 2). Examination of 
residual and Q-Q plots reveals reasonable fits for both the 
binomial and positive models (Suppl. Fig. 2) (online only). 


White shrimp 


White shrimp occurred in 51% of tows (Table 1) and were 
the most abundant penaeid shrimp species caught in CTS 
sampling. White shrimp were most abundant during the 
fall, relative to the other seasons (Fig. 4A). White shrimp 
were spatially distributed throughout the region but were 
most abundant in waters of northern Florida, with smaller 
peaks in abundance off South Carolina and New River, 
North Carolina (Fig. 4B). White shrimp were most abun- 
dant in bottom temperatures ranging between 20.1°C and 
26.5°C (Fig. 4C) and in waters with bottom salinities less 
than 34.9 (Fig. 4D). 
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Nominal abundance of white shrimp (Litopenaeus setiferus), by (A) season, (B) stratum, (C) bottom temperature, and 
(D) bottom salinity, based on data from sampling off the Atlantic coast of the southeastern United States from 1990 
through 2019, with design-based imputation of missing values. (E) Nominal and standardized relative abundances of 
white shrimp estimated by using the design-based method and delta-generalized linear models, respectively, are pro- 
vided by year. Error bars indicate the standard errors of the means. The horizontal gray line indicates the long-term 
mean of 1 for relative abundance over the full period of the time series. 


Nominal abundance of white shrimp was generally 
below, or relatively similar to, the mean for the time series 
from 1990 through 2003 (Fig. 4E). From 2004 through 
2019, nominal abundance of white shrimp was variable 
but with a generally increasing trend compared with the 
trend of the time series mean. Variability in nominal 
abundance generally increased as abundance increased. 

The trend in delta-GLM standardized abundance of 
white shrimp was similar to that in nominal abundance 
but with lower variability (Fig. 4E, Suppl. Table [online 
only]). Although the error bars overlap, the only notable 
difference between nominal and standardized estimates 
was in the abundance estimate for 2019, when the stan- 
dardized abundance of white shrimp was higher than the 
nominal abundance. All covariates were statistically sig- 
nificant in predicting the presence (binomial model) and 
abundance (positive model) of white shrimp (Table 2). 
Examination of residual and Q-Q plots reveals reason- 
able fits for both delta-GLMs (Suppl. Fig. 3) (online only). 


All species 


For the 3 species examined in this study, the standardiza- 
tion of abundance through the use of delta-GLMs primar- 
ily affected the variability in annual estimates, reducing 
the average variability in nominal values for Atlantic 
croaker by 3%, for bluefish by 90%, and for white shrimp 
by 75% (Figs. 2E, 3E, and 4E). The magnitude of the vari- 
ability of abundance estimates generally changed with the 
extent of spatial and temporal distribution and percent 
positive rate (i.e., the percentage calculated as the num- 
ber of tows in which a given species was collected divided 
by the total number of tows conducted per year) among 
the 3 species examined. The Atlantic croaker had the wid- 
est distribution and highest occurrence of the 3 species, 
with the smallest annual nominal errors. The bluefish had 
the narrowest distribution and lowest occurrence of the 
3 species and had an average normalized nominal stan- 
dard error over 2 times higher than that for the Atlantic 
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croaker (0.60 for bluefish and 0.26 for Atlantic croaker) and 
just slightly higher than that for the white shrimp (0.55). 
Standardization had little effect on annual abundance 
estimates, with the exception of increasing abundance esti- 
mates in 2019 for Atlantic croaker and in 2001 for bluefish 
(based on error bars that do not overlap; Figs. 2E and 3E). 
The years with strata that were missed in sampling (1990, 
2013, 2018, and 2019) did not have noticeable differences 
between nominal and standardized abundance estimates, 
with the exception of 2019 for Atlantic croaker and white 
shrimp (although the latter had ovenavuIte error bars) 
(Table 1, Figs. 2E and 4E). 

The effect of scenarios of eit survey modifications 
on standardized abundance was not consistent among 
the 3 species examined in this study (Fig. 5, A-C). The 
Atlantic croaker had the most variable change to index 
values of the 3 species, with the scenario in which the 
dropped season was systematically rotated generally 
tracking most closely with the abundance index from the 
delta-GLM analysis that used the current survey design, 
which includes all seasons (Fig. 5A). The effect of each 
of the 3 scenarios in which a particular season was con- 
sistently dropped on the abundance of Atlantic croaker 
depended on the year examined, likely because the annual 
abundance is informed by the mean response to the other 
seasons, which itself is inconsistent among years, seasons, 
and strata. The effect of randomly selecting a season to be 
dropped was also relatively variable for Atlantic croaker 
and resulted in noticeable deviations from the values in 
the abundance index produced in the model analysis that 
included all seasons. 

All 5 scenarios had minimal effects on the standardized 
abundance index for both bluefish and white shrimp rela- 
tive to the model for all seasons (Fig. 5, B and C), and this 
outcome is most likely related to the more restricted (and 
predictable) spatiotemporal distribution of abundance for 
these 2 species relative to that of the Atlantic croaker. The 
few examples of noticeable deviations from the full index 
occurred for bluefish and white shrimp under the scenar- 
ios in which a single season was consistently not sampled 
and for white shrimp under the scenario in which the 
dropped season was rotated. The variability in the scenar- 
ios for survey design modifications differed among species 
and years, but not in a consistent way (Table 3). 


Discussion 


As a long-term fishery-independent survey, the CTS pro- 
vides abundance and life history data for a variety of spe- 
cies of management interest. A stratified-random sampling 
design is used to conduct the CTS, with survey stations 
allocated to minimize variability in catch abundance esti- 
mates. Historically, annual sampling at the vast majority 
of allocated stations was completed and sampling occurred 
in all strata. However, in recent years (i.e., 2018 and 2019), 
allocated stations and full strata were left unsampled 
because of the use of an aging vessel, mechanical fail- 
ures, stagnant funding and increased costs, and prolonged 


periods of weather conditions above levels deemed safe for 
sampling. In particular, strata at the northernmost and 
southernmost portions of the survey area were those most 
likely to remain unsampled. 

In this study, we investigated whether or not our 
inability to complete the full spatial scope of the survey 
in each season affected the accuracy of abundance esti- 
mates by comparing nominal (design-based) and stan- 
dardized (delta-GLM) estimates of abundance for Atlantic 
croaker, bluefish, and white shrimp. Surprisingly, only one 
major difference between estimates produced with these 
2 methods was observed: one for Atlantic croaker for a 
year with incomplete sampling. Other than this instance, 
there were almost no discernible differences in estimates 
of abundance for the years with incomplete sampling rel- 
ative to estimates for years with complete sampling. One 
reason for the lack of differences in estimates may be the 
use of average long-term values for nominal abundances 
for years with incomplete sampling versus the use of the 
standardization technique in which average values are 
applied within each model cell to account for missing 
data. Although standardization may involve the use of 
an improved average estimate because of the inclusion 
of more covariates than the relatively simple imputation 
method used to produce the nominal estimate (Walters, 
2003), any loss of precision in an index resulting from 
decreased or lost effort can be balanced if stratification is 
incorporated into the design of a given survey and because 
of the similarity of using averages as empirical informa- 
tion (Xu et al., 2015). 

In the case of the CTS, the stratified-random sampling 
design has been consistent over time, allowing minimiza- 
tion of loss of information when sampling was incomplete. 
This assumption that the design will balance incomplete 
sampling, however, may not remain true if completion of 
sampling continues to be eroded. It is currently unknown 
what the effect of incomplete sampling may be for species 
that are less common or less ubiquitous than the 3 species 
examined in this study. 

The relationships between environmental covariates 
included in the standardization models and species abun- 
dance varied among the species examined in this study. The 
association between abundance and bottom temperature 
varied among species. There was little difference in abun- 
dance among temperatures for Atlantic croaker (Fig. 2C), 
but bluefish were most associated with temperatures 
<20.1°C, likely related to their presence in the northern 
strata of the study area in the spring and fall (Fig. 3, B 
and C). Abundance of white shrimp was less influenced by 
temperature than estimates for bluefish, but abundances 
generally increased with increasing temperatures for white 
shrimp, most likely reflecting the warmer waters in south- 
ern strata where white shrimp were abundant (Fig. 4, B 
and C). Salinity was less influential on species-specific 
abundance: there was no significant difference in abun- 
dance of bluefish among salinity levels, Atlantic croaker 
were significantly more abundant in waters with higher 
salinities (=34.9) than in those with lower salinities (<34.9), 
and abundance of white shrimp followed the opposite trend. 
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Figure 5 


Abundance estimates for (A) Atlantic croaker (Micropogonias undulatus), (B) blue- 
fish (Pomatomus saltatrix), and (C) white shrimp (Litopenaeus setiferus) stan- 
dardized with generalized linear models for the 5 scenarios of potential design 
modifications of the Coastal Trawl Survey of the Southeast Area Monitoring and 
Assessment Program—South Atlantic: dropping spring (April and May), summer 
(July and August), or fall (September through November), rotating which season is 
dropped, or randomly dropping a season. The full time series was used in models 
with the survey design modifications applied to sampling conducted between 2014 
and 2019 off the Atlantic coast of the southeastern United States. The black line in 
each panel indicates the standardized estimates of relative abundance produced 
under the current survey design, in which no season is dropped, with error bars for 
standard errors of the annual means. The dashed orange line in each panel indi- 
cates the long-term mean of 1 for relative abundance over the full period of the time 
series (1990-2019). 
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Table 3 


Standard errors of standardized abundance estimates for (A) Atlantic croaker (Microp- 
ogonias undulatus), (B) bluefish (Pomatomus saltatrix), and (C) white shrimp (Litope- 
naeus setiferus) off the Atlantic coast of the southeastern United States between 2014 
and 2019 for the 5 scenarios of potential design modifications of the Coastal Trawl Sur- 
vey of the Southeast Area Monitoring and Assessment Program—South Atlantic: drop- 
ping spring (April and May), summer (July and August), or fall (September—November), 
randomly dropping a season, and rotating which season is dropped. For comparison, 
standard errors of abundance estimates produced under the current survey design, 
with no season dropped, are provided. 


All 
Year seasons Spring Summer Fall Random Rotated 


Atlantic croaker 
2014 58.54 109.40 27.66 
2015 63.67 46.00 39.89 
2016 40.72 29.07 34.23 
2017 34.36 34.77 27.46 
2018 127-0) 50.56 129.87 
2019 242.09 176.95 82.46 
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Bluefish 
2014 
2015 
2016 
2017 
2018 
2019 

White shrimp 
2014 
2015 
2016 
2017 
2018 
2019 


Population levels for Atlantic croaker are driven by 
the combination of temperature and salinity, with winter 
temperatures being primarily influential on survivorship 
and salinities above 10 potentially mitigating the effects 
of extremely low temperatures on survivorship (Lank- 
ford and Targett, 2001; Hare and Able, 2007). Bluefish are 
known to overwinter between North Carolina and Florida, 
and they generally migrate northward in the spring as 
water temperatures increase and then migrate southward 
in the fall as water temperatures decrease (Wilk’; Fahay 
et al., 1999). In previous studies, the low temperature tol- 
erance of juvenile bluefish was found to be 13—15°C (Hare 
and Cowen, 1996), and their growth slowed at water tem- 
peratures greater than 24°C (Hartman and Brandt, 1995). 
It is possible that the binary binning of bottom salinity 
meant the relationship could not be resolved for bluefish 
or that bluefish were present at times with more variable 
salinities (e.g., variable run-off in fall due to seasonal 


’ Wilk, S. J. 1977. Biological and fisheries data on bluefish, Poma- 
tomus saltatrix (Linnaeus). Natl. Mar. Fish. Serv., Northeast 
Fish. Sci. Cent., Sandy Hook Lab. Tech. Ser. Rep. 11, 44 p. [Avail- 
able from website.] 


0.11 0.14 0.11 
0.15 0.18 0.19 
0.09 0.08 0.08 
0.05 0.05 0.04 
0.12 0.07 0.12 
0.22 0.14 0.18 


1.11 2.98 1.78 
4.54 4.00 4.23 


10.74 6.80 7.08 
11.58 10.82 11.42 


4.95 3.62 4,76 


10.19 9.30 9.10 


rains and tropical systems). The reduced abundance of 
white shrimp above a salinity of 34.9 likely reflects their 
high abundances around large estuarine systems with 
relatively high freshwater or estuarine outflow. Results 
from previous studies support the notions that abundance 
of adult white shrimp is positively correlated with water 
temperature and negatively correlated with salinity and 
that mortality occurs at temperatures below 10°C (Joyce, 
1965; SAFMC, 1981; Diop et al., 2007; Fowler et al., 2018). 

In general, standardization decreased variability relative 
to design-based estimation of abundance in 18 of 30 survey 
years (60%) for Atlantic croaker, in 12 survey years (40%) 
for bluefish, and in 16 survey years (53%) for white shrimp 
(Figs. 2E, 3E, and 4E). For all 3 species examined in this 
study, the primary variables of season and stratum in the 
design-based GLMs accounted for more of the deviance 
explained by GLMs than either bottom temperature or salin- 
ity and often explained levels of deviance similar to those 
explained by year (Table 2). Both the spatial and temporal 
elements of the survey design enabled us to find that those 
species with more restricted distributions have greater vari- 
ability in estimates of abundance. This pattern of small dis- 
tribution and larger variability may be partially accounted 
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for through the nominal estimation used in this study but 
is likely more adequately addressed by using the standard- 
ization method that can incorporate other drivers such as 
environmental conditions (Figs. 2E, 3E, and 4E). 

The bluefish had the most uneven spatial distribution 
of the 3 species examined and had some of the highest 
annual standard error estimates, but the white shrimp 
was the species most unevenly distributed among seasons 
and had moderate standard error estimates (Figs. 3B, 3E, 
4B, and 4E). Carruthers et al. (2011) suggested that the 
choice to impute missing data rather than to standard- 
ize by using conventional GLM techniques should be 
dependent on the spatial dynamics of a population and 
the spatial grain of the missing data. This choice between 
methods is especially important if there are area and time 
interactions and the unobserved areas do not adhere to 
the mean abundance trend. Use of imputation tends to 
lead to greater bias than use of a GLM, unless the spatial 
scale of imputation is relatively fine and the imputed data 
are representative of the overall population. 

In our study, we chose to impute a mean based on the full 
time series, although imputing a mean based on data from 
only recent years may have been more representative of 
more recent trends for each species. However, for bluefish, 
using stratum in the imputation method is more appropri- 
ate than using larger latitudinal bins because stratum has 
the finest spatial scale consistent with the survey design 
(station cannot be used as a variable because not every 
station is sampled each year). Because of the uneven spa- 
tial distribution of bluefish, future work should determine 
if the use of station for imputation is more appropriate 
than the use of stratum. One caveat to using station is 
that, because not every station is selected for sampling 
each year, the use of station may result in a decrease in 
sample sizes that inform calculations. 

There were only 2 observed differences in abundance 
estimates between the nominal and standardization meth- 
ods as revealed by error bars not overlapping for Atlantic 
croaker in 2019 and for bluefish in 2001. The year 2019 
was the only year of incomplete sampling during which 
the 3 northernmost strata were not sampled in both spring 
and fall. In general, these strata are also the ones that typi- 
cally have moderate to high abundances of Atlantic croaker 
(Fig. 2B). Not surprisingly, for 2019, the nominal abundance 
estimate was markedly lower than the standardized value. 
The difference between the nominal and standardized 
abundance estimates is likely a result of the standardiza- 
tion method better accounting for the effects of covariates 
on abundance, for the general increasing trend in 2019, 
and for the typically higher catches in these strata than 
the design-based approach, even with the average nomi- 
nal abundance imputed as a proxy for the missing data. By 
comparison, there was no change in the estimates for either 
bluefish or white shrimp in 2019. Catches of white shrimp 
in survey tows in the northernmost strata were generally 
low (Fig. 4B); therefore, missing data for these strata likely 
did not have a great effect on the annual abundance esti- 
mate from either calculation. Generally, the highest abun- 
dances observed for bluefish occurred in the 3 northernmost 


strata and in spring and fall (Fig. 2, A and B). This consis- 
tent positive relationship between stratum (latitude) and 
abundance likely allowed those consistently high catches 
to be captured in the calculations of both the nominal and 
standardized methods, resulting in very similar estimates 
for 2019. In contrast to survey effort in 2019, all allocated 
stations were sampled in 2001. 

Two of the few exceptions to improvement in accuracy 
(i.e., reduced error) of the standardized abundance esti- 
mates were for Atlantic croaker in 2019 and for bluefish 
in 2001. For Atlantic croaker in 2019, the estimate of 
abundance increased significantly following standard- 
ization, and the error rate increased almost 3-fold. Strata 
were missed in sampling or not sampled completely in 
both spring and fall in 2019, the most variable and least 
abundant seasons in the full time series for Atlantic 
croaker, respectively. Abundance of Atlantic croaker in 
fall 2019 was higher than in any other previous fall but 
also was more variable, likely a result of not sampling in 
the northernmost strata where this species is abundant 
but the abundance is variable and of not sampling in the 
southernmost strata where this species is least abundant, 
reducing the number of strata in which there were zero 
or near zero encounters. In addition, the use of imputa- 
tion in the design-based calculation did not account for 
uncertainty in the imputed estimates for multiple seasons 
and strata and, therefore, reduced variability for 2019 by 
default (Carruthers et al., 2010). In contrast, errors for all 
seasons and strata in 2019 continued to be incorporated 
when the standardization technique was used, potentially 
leading to the higher variability in estimates for that year. 
Further comparisons should examine variability of esti- 
mates based on imputed data. . 

Regarding bluefish in 2001, there is no record of 
unusual circumstances in CTS historic data for that year; 
therefore, there may have been unusual conditions that 
were not accounted for in the current analyses. If so, these 
conditions were likely specific to bluefish because the esti- 
mates for the other 2 species in 2001 were similar and had 
decreased standardized variability as expected. 

Results of the exercise of retrospectively dropping sea- 
sons for the years 2014-2019 indicate that the likely 
effects of major modifications to the survey would be 
species-specific (i.e., not all species abundance estimates 
remained consistent among all modification scenarios). 
Removing whole seasons did not change the pattern of 
increased abundance of Atlantic croaker in 2019 for most 
scenarios: dropping fall, dropping summer, and systemat- 
ically rotating which season was dropped from analysis 
(Fig. 5A). By comparison, removing spring or randomly 
dropping seasons between 2014 and 2019, decreased the 
standardized abundance estimate for Atlantic croaker in 
2019 compared with the standardized estimate produced 
with data for all seasons. Interestingly, dropping seasons 
did not change the pattern of abundance for bluefish or 
white shrimp (Fig. 5, B and C). Both species had lower 
abundance estimates and more variability compared with 
estimates for Atlantic croaker, but with standardization 
we were able to successfully predict abundances for these 
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species even with the loss of full seasons. Therefore, there 
may be an interaction between annual abundance of each 
species, species distributions or predictability, strata or 
seasons that are sampled, and effects of standardization 
on variability and estimates of abundance that need to be 
considered before making changes to the survey design. 

Fishery-independent indices of abundance generally are 
essential for most stock assessment models currently in 
use, are preferred over fishery-dependent indices, and are 
often upweighted relative to other indices because of their 
adherence to standardized designs and resistance to tem- 
poral changes in catchability or sampling systems (Wilberg 
et al., 2009). For the CTS, delta-GLM indices of abun- 
dance indicate that standardization decreased variability 
in annual abundance estimates for all species examined 
and potentially “corrected” data from the year in which 
sampling was most incomplete for one of the species most 
commonly encountered in the survey, the Atlantic croaker. 
On the basis of the improvements in variability alone, we 
recommend moving forward with standardization for spe- 
cies for which long-term time series are available. Because 
of the advanced age of the vessel used for surveys and the 
trend of more extreme weather events, standardization 
may become even more important in future years (and for 
more species) as completion of CTS sampling is likely to 
continue to be challenging. 

Given that the survey design needs to be modified to 
account for these potential constraints on sampling com- 
pletion, that the effect on variability among modification 
scenarios was not informative, and that the best perform- 
ing model varied among species and years (Table 3), we 
also recommend rotating the season that is not sampled, 
rather than any of the other modification scenarios exam- 
ined in this study, if modifications become necessary. The 
effect of a given season on the abundance index is species 
specific, and results of this study indicate that rotation of 
the dropped season is the best approach for a multispecies 
survey if sampling effort needs to be reduced. 


Conclusions 


In this study, as expected, a common standardization tech- 
nique (the delta-GLM method) improved estimation of vari- 
ability in abundance time series from a fishery-independent 
survey over that of the simpler approach of imputation of 
average abundance values for Atlantic croaker, bluefish, 
and white shrimp. These improvements in uncertainty 
become increasingly important as costs of conducting sur- 
veys increase and funding for long-term surveys stagnates, 
as has been experienced for the CTS. The results of this 
study indicate that survey programs that cannot complete 
sampling as designed in a given year could deal with that 
incompletion through techniques that involve a delta-GLM. 
We also found that, when funding is severely limited, 
rotating the season that is not sampled may be the most 
appropriate approach for multispecies surveys. For surveys 
without a seasonal structure, incorporating a time element, 
such as month, a posteriori and rotating which months are 


prioritized may be a way to cope with limited funding and 
sampling effort as well, although it would need to be exam- 
ined on a case-by-case basis. 

In addition, for the CTS, over 24 species currently are of 
management interest and are sampled sufficiently to pro- 
duce indices of abundance. For this high variety of species, 
more work needs to be done to examine the effects of sur- 
vey design modifications because the results of this study 
were not consistent among species. The delta-GLM method 
is just one of many standardization techniques that can be 
employed. Other types of models, such as those that use con- 
tinuous covariates or other error structures, may be more 
appropriate or better at improving variability estimation 
than the analyses for which results are presented herein. 


Resumen 


Los indices de abundancia procedentes de estudios inde- 
pendientes de la pesca son los preferidos en las evaluaciones 
de poblaciones debido a sus sélidos disefios cientificos que 
minimizan la incertidumbre y el sesgo. Cuando el muestreo 
no se ajusta al disefio, los investigadores emplean técnicas 
como la imputacién o la estandarizacién para mejorar la 
precisién y reducir el sesgo. Examinamos 2 métodos para 
ajustar el muestreo incompleto en la Prospeccién de Arra- 
stre Costero (CTS) del Programa de Monitoreo y Evalu- 
acion del Area sureste del Atlantico Sur para 3 especies 
comunmente encontradas en el muestreo de la prospec- 
cion, la corvina del Atlantico (Micropogonias undulatus), la 
anjova (Pomatomus saltatrix) y el camar6én blanco (Litope- 
naeus setiferus): la imputacion de datos faltantes con base 
en el disefio y la estandarizacion a través del enfoque del 
modelo lineal generalizado delta. Ademas, determinamos 
el efecto de modificar el componente estacional del disefio 
de la prospeccién mediante una simulacion retrospectiva. 
Para las 3 especies, la estandarizacio6n mejor6 la precisién 
de las estimaciones de abundancia anual con relaci6n a 
los valores estimados con el método basado en el disefio. 
Cuando un estrato omitido en el muestreo coincidia con 
un area o época de alta variabilidad para una especie 
(p. ej., 2019), la estandarizacién no mejoré6 la precisién con 
respecto al método basado en el disefio. Los resultados 
del andlisis de los efectos de la exclusién de temporadas 
enteras, debido a problemas de financiacién o logisti- 
cos, indican que la rotacién de la temporada a excluirse 
fue el mejor enfoque para equilibrar las caracteristicas 
de cada especie. En general, recomendamos el enfoque de 
estandarizacién para contabilizar los datos faltantes en 
las series temporales de la CTS. 
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Abstract—Heavy rain can decrease 
salinity and increase turbidity of the 
water in coastal areas, negatively affect- 
ing the development of organisms, par- 
ticularly during their early life stages. 
In this study, the effects of salinity and 
turbidity on embryos of the bamboo sole 
(Heteromycteris japonicus) were evalu- 
ated to improve understanding of its tol- 
erance to global climate change. Three 
experiments were carried out over a 7-d 
period. In the first experiment, embryos 
of bamboo sole were exposed for 3 h to 1 
of 6 salinity levels (14-34). Low salinity 
levels (14 and 18) resulted in signifi- 
cantly shorter total lengths of newly 
hatched larvae in comparison with 
larval sizes in treatments with higher 
salinities, but no significant differences 
were observed in hatching rate and lar- 
val survival rate among treatments. In 
the second experiment, embryos were 
exposed to turbidities of 0, 100, 300, 
500, and 700 nephelometric turbidity 
units for 3 h. Turbidity significantly 
decreased hatching rate, survival rate, 
and total length and increased onset 
hatching time and percentage of abnor- 
mality. In the third experiment, embryos 
were exposed to different combinations 
of salinity and turbidity. The interaction 
effect of salinity and turbidity on total 
length of newly hatched larvae was sig- 
nificant. These findings indicate that 
embryo development of bamboo sole 
was more affected by changes in turbid- 
ity than by changes in salinity. 
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Flatfish species are found globally in 
marine habitats and are important 
food fish to humans (Gibson et al., 
2015). However, the global catch of 
flatfish in fisheries has decreased grad- 
ually since the 1970s, from about 3 mil- 
lion metric tons (t) to below 2 million t 
in the 2000s, and is predicted to decline 
by 20% by the year 2100 (Cheung and 
Oyinlola, 2018). 

Although overfishing has resulted 
in decreases in flatfish populations, 
there is another related reason for these 
reductions: climate change (Cheung 
and Oyinlola, 2018). In fact, in many 
studies, fluctuations in fish populations 
have been found to be closely associated 
with climate change (Lehodey et al., 
2006; Brander, 2010; Holsman et al., 
2012; Barange et al., 2014). Baudron 
et al. (2014) also reported that negative 
changes in environmental conditions 
of habitats due to climate change can 
directly affect fish population dynam- 
ics through physiological effects. 

One of the effects of climate change 
in coastal areas is the increased inten- 
sity and frequency of heavy rainfall, 
and those changes in rainfall reduce 


salinity and increase turbidity as a 
result of runoff from the land and riv- 
ers (Thrush et al., 2004; Poloczanska 
et al., 2009). According to Allen and 
Pechenik (2010), salinity in the upper 
water column may fluctuate from 30 to 
15 in short periods during heavy rain- 
fall events, leading to pelagic fish eggs 
likely being exposed to low salinity 
during development. In a recent study, 
Fukuda et al. (2021) found that, after 
heavy rains, the salinity level was less 
than 28 in Hakata Bay, which is located 
in the city of Fukuoka on the northern 
end of Kyushu in Japan. In the natural 
environment, turbidity is usually less 
than 50 nephelometric turbidity units 
(NTU) (Boyd, 2015); however, heavy 
rainfall events cause high runoff to 
flow more directly into coastal areas, 
resulting in turbidity values that can 
be greater than 50 NTU. During heavy 
rains, the turbidity of outflow can 
exceed 130 NTU (Zhou et al., 2015) 
and even reach levels as high as 467 
NTU (Zheng et al., 2016), but it rarely 
exceeds 1000 NTU (Boyd, 2015). 
Seawater salinity and turbidity are 
2 important ecological variables that 
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have individual or interactive effects on development, 
growth, and survival of fish at any life history stage (Boeuf 
and Payan, 2001; Hasenbein et al., 2013; Phan et al., 
2019). In particular, drastic changes in coastal salinity 
and turbidity, such as those due to extreme weather 
events, can profoundly affect egg fertilization and incu- 
bation, yolk sac absorption, early embryogenesis, and lar- 
val survival (Berlinsky et al., 2004; Bilotta and Brazier, 
2008). For example, the short-term effect of low salinity 
(<3) was to lower embryonic survival and egg hatching in 
a study of Pacific herring (Clupea pallasii) (Mikhailenko, 
2000). Gray et al. (2012) found that short-term experi- 
mental elevation of turbidity had a significant effect on 
the hatching success rate of embryos of the spotted gar 
(Lepisosteus oculatus), a species that was then consid- 
ered threatened and is now considered endangered in 
Canada (COSEWIC, 2015). Phan et al. (2019) reported 
that increased larval abnormality was recorded in bas- 
tard halibut (Paralichthys olivaceus) when embryos were 
exposed to high turbidity (>100 NTU) or low salinity 
(<30). Although prior to the study we describe herein the 
short-term effects of salinity and turbidity on the embry- 
onic stage of bamboo sole were unknown, these findings 
indicate how the bamboo sole (Heteromycteris japonicus) 
might respond to climate change. 

The bamboo sole is a small soleid flatfish distributed 
along southern Japan, in the Yellow Sea, and in the East 
China Sea (Nakabo, 2002). It is a common species found 
on sandy seabeds in the Seto Inland Sea of Japan (Yama- 
moto and Katayama, 2013). Although the bamboo sole is 
a noncommercial species, it plays an important ecolog- 
ical role in coastal areas. As reported by Jawad (2021), 
the early life stages of many noncommercial species are 
essential components of the diets of important commer- 
cial species. The spawning season for bamboo sole is 
from May through June, and their pelagic eggs are often 
caught in estuaries or bays (Yamamoto and Katayama, 
2013), where salinity and turbidity fluctuate because 
of heavy rainfall events. Therefore, the objective of this 
study was to evaluate the effects of salinity and turbidity 
changes on various characteristics of the early stage of 
bamboo sole. Specifically, we examined the hatching time 
and hatching rate of embryos, the percentage of abnor- 
mality, size, and yolk sac volume of newly hatched larvae, 
and the survival rate of larvae. 


Materials and methods 
Collection and preparation of fertilized eggs 


Embryos of bamboo sole were obtained in July 2019 from 
the National Research Institute of Aquaculture, which 
is part of the Fisheries Research Agency of Japan and 
located in the Mie Prefecture. The fertilized eggs were 
collected and transferred to the Shallow Sea Aquaculture 
laboratory, Mie University, in oxygen-filled plastic bags. 
Then, eggs were placed in a 5-L glass tank at a salinity of 
34 and a temperature of 23°C, with gentle aeration. Only 


viable (buoyant) embryos were used for our experiments. 
The fertilization rate of eggs was 80.2%, and the mean 
diameter of eggs was 1.02 mm (standard error of the 
mean [SE] 0.04). For experiments, we created seawater 
at different levels of salinity by diluting artificial seawa- 
ter powder (LIVESea Salt’, Delphis Co. Ltd., Itami City, 
Japan) with distilled water and checking the salinity with 
a handheld refractometer (Master-S/Mill Alpha, Atago 
Co. Ltd, Tokyo, Japan). To achieve the desired turbidity, 
kaolin clay (with the chemical composition Al,Si,0;(OH), 
and a particle diameter of approximately 0.4 pm) was 
dissolved in water with a salinity of 34. Turbidity was 
measured by using a HACH DR/850 portable colorimeter 
(Hach Co., Loveland, CO). 


Experimental design 


The experimental conditions were set up to mimic, for a 
short period, the stressful conditions that organisms typ- 
ically are exposed to in bays, estuaries, and coastal areas 
during heavy rains, and then conditions were returned to 
normal (Phan et al., 2020). Salinity and turbidity are the 
2 major variables for water quality in coastal areas of the 
ocean during and shortly after heavy rainfall events (Phan 
et al., 2020). In this study, 3 experiments were conducted 
to examine the following: 1) effect of low salinity levels (6 
treatments with salinities of 34 [control], 30, 26, 22, 18, 
and 14) on embryos of bamboo sole, 2) effect of high tur- 
bidity levels (5 treatments with turbidities of 0 [control], 
100, 300, 500, and 700 NTU) on embryos of bamboo sole, 
and 3) combined effect of salinity and turbidity (8 treat- 
ments with combinations of salinities of 14 and 18 and 
turbidities of 100, 300, 500, and 700 NTU) on embryos of 
bamboo sole. 

All experiments were conducted in 6-well plastic micro- 
plates, so that each treatment had 6 replicates. When 
embryos developed to the gastrula stage, 13 normal 
embryos were removed from the 5-L tank and randomly 
allocated into each of the 6 wells of a plastic microplate, 
with each well containing 5 mL of artificial seawater (with 
a salinity of 34). The embryos were exposed to single and 
combined salinity and turbidity levels for 3h. Then embryos 
were transferred to 100-mL glass beakers with seawater 
at a temperature within 1°C of the target temperature of 
23°C, a salinity of 34, and a turbidity of 0 NTU under a 
natural photoperiod (embryos were exposed to 12 h of day- 
light per day) in an electric incubator. Embryos remained 
in the beakers in the incubator until the end of the exper- 
iment. Gradual changes in salinity or in combined condi- 
tions of salinity and turbidity were achieved by using an 
automatic pipette and a water exchange rate of 1 mL until 
the desired environmental condition was reached. During 
the incubation period, 50% of the seawater in all treat- 
ments was exchanged daily, and dead embryos and larvae 
were removed by using a wide-mouthed pipette. Embryos 
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were determined to be dead if morphology had collapsed 
(embryos appeared white or opaque) or the heartbeat had 
stopped. Seawater in the beakers had dissolved oxygen of 
5.2-5.7 mg/L and pH of 8.0—8.5. Larvae at 3 d post-hatch 
were fed with live rotifers, specifically with Brachionus 
plicatilis (4 rotifers/mL). All experiments were carried out 
for a period of 7 d. 


Measurements and calculations 


Onset hatching time, hatching rate, percentage of abnor- 
mality, survival rate, total length, and yolk sac volume of 
embryos and larvae were recorded. Onset hatching time 
was determined as the total number of hours from fertil- 
ization to the time of hatch. Hatching rate was calculated 
on the basis of the number of hatching embryos compared 
to the number of stocked fertilized eggs. The larval survival 
rate was determined as the total number of live larvae at 
5 d post-hatch. The percentage of abnormality at 24 h post- 
hatch was calculated by using the numbers of normal and 
abnormal larvae. Individuals were considered abnormal if 
the notochord was bent or the tail tip was deformed. Only 
a few abnormal larvae appeared across all treatments in 
the first experiment; however, all abnormal larvae died by 
24 post-hatch. Therefore, the percentage of abnormality 
was not recorded. Shortly after hatching, 18 larvae from 
each treatment were randomly collected and preserved in 
4% formalin solution, and total length and yolk sac volume 
of newly hatched larvae were measured under an inverted 
microscope (CKX53, Olympus Corp., Tokyo, Japan) with 
an ocular micrometer that had a precision of 10 ym. Yolk 
sac volume (in cubic millimeters) of newly hatched larvae 
was calculated by using this formula: 


Yolk sac volume = - x LH?, 


where L = the length of the yolk sac (in millimeters); and 
H = the height of the yolk sac (in millimeters). 


Statistical analysis 


The software SPSS Statistics (vers. 22.0, IBM, Armonk, 
NY) was used for statistical analyses. All data are 
expressed as means with their SEs. Results were analyzed 
by using 1-way (first and second experiments) and 2-way 
(third experiment) analysis of variance (ANOVA). Before 
analysis, percentages in data were arcsine transformed. 
Significant treatment effects were detected by using the 
Tukey’s honestly significant difference test, a post hoc 
multiple-range test. The significance level was 0.05. 


Results 
Effect of low salinity levels on embryos 


In the first experiment, 468 embryos of bamboo sole were 
exposed to a control level of salinity or to 1 of 5 different 
salinities (78 individuals for each of 6 treatments). Embryos 
had great tolerance to a gradual decrease in salinity. In fact, 
salinity had no effect on onset hatching time (range: 39.3— 
39.7 h) (1-way ANOVA: F[5, 30]=46.39, P>0.05) or on hatch- 
ing rate (range: 91.0-97.4%) (1-way ANOVA: F|5, 30]=0.649, 
P>0.05). Larval total length was positively related to salin- 
ity, ranging from 1.69 mm at a salinity of 34 to 1.61 mm at 
a salinity of 14. The total lengths of newly hatched larvae 
in the treatments with salinities of 14 and 18 were signifi- 
cantly smaller than those of larvae in the control treatment 
with a salinity of 34, but no significant differences in the 
total lengths of larvae were recorded between the treat- 
ments with salinities of 14 and 18 or among treatments 
with salinities of 22, 26, 30, and 34. Yolk sac volume ranged 
from 0.100 to 0.110 mm? across all treatments, but no sig- 
nificant differences were observed among treatments. 

On day 5 after hatching, there were no significant dif- 
ferences in survival rates among all salinity treatments 
(F[5, 30]=1.626, P>0.05) (Table 1). 


Table 1 


Mean measurements of characteristics of the early stage development of bamboo sole (Hetero- 
mycteris japonicus) in laboratory experiments conducted at Mie University in Japan in July 
2019 to examine the effect of salinity on development. Embryos were exposed to 1 of 6 levels 
of salinity for 3 h. After treatment, embryos were transferred to a controlled environment 
and monitored for 7 d. The survival rate of larvae was measured 5 d post-hatch. The different 
letters after values indicate that the means are significantly different (P<0.05) among treat- 
ments. Standard errors of the means are provided in parentheses. 


Yolk sac 
volume (mm?) 


Hatching 
rate (%) 


Survival Total length 
rate (%) (mm) 


Onset hatching 
time (h) 


Salinity 
treatment 


97.4 (1.6)? 
93.6 (3.1)? 
92.3 (4.0)? 
93.6 (2.4)° 
91.0 (2.4)? 
01.0137)" 


80.0 (2.6)* 
78.3 (4.0)* 
75.0 (2.2)* 
TREY 
73.3 (3.3)* 
70.0 (2.6)* 


1.69 (0.01)? 
1.68 (0.01) 
1.67 (0.01)* 
1.67 (0.01)? 
1.62 (0.01)? 
1.61 (0.01)? 


0.110 (0.004)* 
0.100 (0.003)* 
0.101 (0.004)* 
0.103 (0.004)* 
0.105 (0.004)* 
0.103 (0.003)* 


34 (control) 39.3 (0.09)? 
30 39.3 (0.07)* 
26 39.5 (0.08)? 
29 39.5 (0.10)? 
18 39.7 (0.18)? 
14 39.5 (0.11)* 
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Effect of high turbidity levels on embryos 


In the second experiment, 390 embryos of bamboo sole 
were exposed to a control level of turbidity or to 1 of 4 dif- 
ferent turbidities (78 individuals for each of 5 treatments). 
Turbidity had a significantly positive effect on the onset 
hatching time of bamboo sole eggs (range: 39.3—41.8 h) 
(1-way ANOVA: F[4, 25]=84.932, P<0.05) (Table 2). The 
times to onset hatching in turbidities of 300 and 500 NTU 
were significantly longer than the time in the control treat- 
ment with a turbidity of 0 NTU and in the treatment with a 
turbidity of 100 NTU and were significantly shorter than the 
time in the treatment with a turbidity of 700 NTU (P<0.05). 
Hatching rate was also affected by turbidity (range: 85.9— 
96.2%) (1-way ANOVA: F[4, 25]=4.667, P<0.05). The lowest 
hatching rate was observed in the treatment with the tur- 
bidity of 700 NTU and was significantly lower than the rate 
recorded in the control treatment and the treatment with a 
turbidity of 100 NTU. However, no significant differences in 
hatching rate were found among treatments with the con- 
trol, 100-, 300-, and 500-NTU turbidities. As for abnormality 
rates, the results of statistical analysis indicate that there 
were significant differences between the control treatment 
and the other treatments (with turbidities of 100, 300, 500, 
and 700 NTU), but there were no significant differences 
among those other treatments. 

The total length of newly hatched larvae significantly 
decreased with increasing turbidity levels. The size of 
larvae in the treatment with a turbidity of 700 NTU was 
significantly lower than the size of larvae in the other treat- 
ments. No significant differences were recorded among 
the treatments with turbidities of 100, 300, and 500 NTU 
(P>0.05). Yolk sac volume in the treatment with 700-NTU 
turbidity was significantly smaller than that in the control 
treatment (P>0.05), but significant effects were not found 
among the control treatment and the treatments with tur- 
bidities of 100, 300, 500 NTU. At the end of the experiment, 
the survival rate of larvae in the treatment with a turbidity 


of 700 NTU was significantly different from the rate in 
the control treatment, but no significant differences were 
found in the larval survival rate between the treatments 
with turbidities of 500 and 700 NTU or among the control 
treatment and the treatments with turbidities of 100, 300, 
and 500 NTU. 


Combined effect of salinity and turbidity on embryos 


In the third experiment, 624 embryos of bamboo sole 
were exposed to 1 of 2 levels of salinity and to 1 of 4 lev- 
els of turbidity (78 individuals for each of 8 treatments). 
Although both salinity and turbidity affect the hatching 
of embryos of bamboo sole, turbidity had a greater effect 
on embryo development and larval survival rate than 
salinity (Tables 3 and 4). The results of 2-way ANOVA 
analysis indicate that turbidity had an effect on onset 
hatching time (F[8, 40]=20.933, P<0.05) but that the 
effects of salinity (F[8, 40]=1.458, P>0.05) and of inter- 
actions of salinity and turbidity (F'[8, 40]=0.972, P>0.05) 
were not statistically significant. The same pattern was 
observed for hatching rates, percentages of abnormality, 
and survival rates. Neither salinity nor the interaction 
between salinity and turbidity affected hatching rate 
(2-way ANOVA: F[8, 40]=0.047, P>0.05), percentage of 
abnormality (2-way ANOVA: F[8, 40]=0.218, P>0.05), 
or survival rate (2-way ANOVA: F[8, 40]=0.278, P>0.05) 
of bamboo sole eggs and larvae. However, turbidity 
had a significant effect on each of these characteristics 
(Table 4). A similar pattern was observed for yolk sac 
volume of newly hatched larvae. Turbidity alone affected 
volume (F[3, 136]=2.768, P>0.05), but there were no sig- 
nificant differences in volume among treatments with 
different salinities (F[3, 136]=0.062, P>0.05) or among 
treatments with interactions of different levels of turbid- 
ity and salinity (F[8, 136]=0.049, P>0.05). Total length 
of newly hatched larvae in all treatments ranged from 
1.53 to 1.65 mm. The results of 2-way ANOVA indicate 


Table 2 


Mean measurements of characteristics of the early stage development of bamboo sole (Heteromycteris japon- 
icus) in laboratory experiments conducted at Mie University in Japan in July 2019 to examine the effect 
of turbidity on development. Embryos were exposed to 1 of 5 levels of turbidity for 3 h. After treatment, 
embryos were transferred to a controlled environment and monitored for 7 d. The survival rate of larvae was 
measured 5 d post-hatch. The different letters after values indicate that the means are significantly different 
(P<0.05) among treatments. Standard errors of the means are provided in parentheses. NTU=nephelometric 
turbidity units. 


Turbidity 
treatment 
(NTU) 


0 (control) 


Onset 
hatching 
time (h) 


39.3 (0.10)? 
39.6 (0.06)? 
40.9 (0.07)° 
40.8 (0.16)? 
41.8 (0.12)° 


Hatching 
rate (%) 


96.2 (1.7) 
94.9 (1.6)? 
91.0 (2.4)” 
88.5 (1.7) 
85.9 (2.4)? 


Percentage of 
abnormality 
(%) 


2.6 (1.6) 
18.8 (3.2) 
19.9 (2.0) 
24.6 (2.5) 
27.0 (3.3) 


Survival 
rate (%) 


S17 Cl,7)" 
7520 (2:2) 
TOOT 
65.0 (3.4) 
56.7 (4.2)° 


Total 
length 


(mm) 


1.69 (0.01)? 
1.64 (0.01)? 
1.63 (0.01)"° 
1.60 (0.01)"° 
1.59 (0.01)° 


Yolk sac 
volume 
(mm°) 


0.110 (0.004)? 
0.094 (0.004)*» 
0.090 (0.006)?” 
0.090 (0.006)?” 
0.089 (0.004)? 


292 


Table 3 
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Mean measurements of characteristics of the early stage development of bamboo sole (Heteromycteris japonicus) 
in laboratory experiments conducted at Mie University in Japan in July 2019 to examine the combined effect of 
salinity and turbidity on development. Embryos were exposed for 3 h to 1 of 8 treatments that combined 2 levels 
of salinity and 4 levels of turbidity. After treatment, larvae were transferred to a controlled environment and 
monitored for 7 d. The survival rate of larvae was measured 5 d post-hatch. Standard errors of the means are 
provided in parentheses. NTU=nephelometric turbidity units. 


Turbidity 
Salinity (NTU) 
18 100 
300 
500 
700 
100 
300 
500 
700 


Onset 
hatching 
time (h) 


40.6 (0.06) 
41.2 (0.12) 
41.4 (0.07) 
41.4 (0.15) 
40.9 (0.12) 
41.3 (0.08) 
41.4 (0.07) 
41.4 (0.05) 


Hatching 
rate (%) 


91.0 (5.0) 
82.1 (2.6) 
80.8 (4.3) 
79.5 (4.3) 
91.0 (3.1) 
83.3 (3.7) 
82.1 (4.7) 
78.2 (3.7) 


Percentage of 
abnormality 


(%) 


15.4 (2.8) 
20.5 (3.2) 
25.6 (3.2) 
26.9 (4.8) 
19.2 (3.3) 
26.9 (3.3) 
26.9 (4.3) 
28.2 (4.3) 


Survival 
rate (%) 


70.0 (7.7) 
58.3 (6.0) 
55.0 (4.3) 
53.3 (4.9) 
71.2 (3.1) 
53.3 (6.7) 
46.7 (6.7) 
48.3 (4.0) 


Total 
length 
(mm) 


1.63 (0.01) 
1.64 (0.01) 
1.65 (0.01) 
1.62 (0.02) 
1.62 (0.02) 
1.61 (0.02) 
1.54 (0.02) 
1.53 (0.02) 


Yolk sac 
volume 
(mm?) 


0.098 (0.004) 
0.095 (0.005) 
0.089 (0.005) 
0.089 (0.005) 
0.099 (0.003) 
0.093 (0.004) 
0.089 (0.003) 
0.087 (0.004) 


significant effects of both salinity (F[3, 136]=18.653, 
P<0.05) and turbidity (F[8, 136]=3.441, P<0.05) on total 
length, with a significant interaction between these vari- 
ables (F'[3, 136]=3.156, P<0.05). 


Discussion 


Shallow coastal areas are some of the most vulnerable sys- 
tems that will be affected by climate change (Phillippart 
et al., 2003). Findings from numerous studies indicate that 
fish in their early life stages were more susceptible than 
adults to fluctuations in environmental conditions and 
that they may be particularly sensitive to climate change 
(Pankhurst and Munday, 2011; Sugisaki and Murakami, 
2018). However, the tolerance levels of marine fish in 
early developmental stages to changes in environmen- 
tal conditions can vary among fish species (Wang et al., 
2002; Shi et al., 2009). Torres et al. (2021) reported that 
organisms exposed to short-term fluctuations in environ- 
mental conditions may experience irreversible damage to 
body structures and eventually suffer mortality at higher 
rates if the range of variation is wider than the range of 
tolerance. In our study, simulating the sudden, short-term 
changes in salinity that occur after a rainstorm, we found 
that exposure to low salinity did not affect onset hatching 
time and hatching rate of eggs and did not affect survival 
rate of larvae after hatching. These findings are similar to 
those of our previous study, in which embryos of Japanese 
seabream (Pagrus major) were subjected to stress from 
low salinity (Phan et al., 2020). 

Salinity tolerance of fish embryos has been reported to 
have improved with advancing development; for example, 
cleaving embryos were more sensitive to salinity change 
than blastulae, and blastulae were more sensitive to 
change in salinity than gastrulae (Farhana et al., 2019). 
Ord (2019) reported that embryos in the late developmental 


stages have more advanced osmoregulatory capacities. In 
the natural environment, coastal fish species are consid- 
ered to be euryhaline. Therefore, their embryos can tol- 
erate a wide range of salinities, as has been found for the 
European seabass (Dicentrarchus labrax) (Conides and 
Glamuzina, 2001), black sea bass (Centropristis striata) 
(Berlinsky et al., 2004), bastard halibut (Phan et al., 2019; 
Manuel et al., 2021), and Japanese seabream (Phan et al., 
2020). 

In our laboratory experiments with bamboo sole, short- 
term exposure to stress from turbidity significantly pro- 
longed onset hatching time and decreased hatching rate 
of embryos and decreased survival rate of larvae after 
hatching. In nature, increased turbidity tends to coincide 
with decreases in dissolved oxygen, which could lead to 
lack of oxygen during development in individuals of many 
fish species (Mueller et al., 2017). In addition, turbidity 
impedes oxygen supply and the normal growth activity of 
embryos as kaolin clay coats egg membranes (Phan et al., 
2020). Del Rio et al. (2019) revealed that oxygen signifi- 
cantly influenced the hatching process in fish embryos, in 
which low levels of oxygen led to delays in embryonic devel- 
opment and in the time of hatching as well as to a decrease 
in hatching success rate. Results of our study indicate that 
the percentage of larval bamboo sole that had abnormal- 
ities was higher when embryos were exposed to turbidity 
levels of 100-700 NTU than when embryos were exposed 
to no turbidity. The morphological abnormalities most 
commonly found in our study were notochord bending and 
tail-tip deformity. An increase in the proportion of larvae 
that were abnormal because of exposure to high turbid- 
ity has also been reported for other fish species, including 
the Pacific herring (Griffin et al., 2009), bastard halibut 
(Phan et al., 2019), climbing perch (Anabas testudineus) 
(Nurulnadia et al., 2020), and Japanese seabream (Phan 
et al., 2020). According to Cahu et al. (2003), such abnor- 
malities are often associated with growth depression and 
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Table 4 


Results from the 2-way analysis of variance in measurements of early 
stage development of bamboo sole (Heteromycteris japonicus) from lab- 
oratory experiments conducted at Mie University in Japan in July 2019 
to examine the individual and combined effects of salinity (S) and tur- 
bidity (T) on development. Embryos were exposed to treatments with 
different levels of salinity and turbidity for 3 h and then monitored for 
7 d. An asterisk (*) indicates that the means are significantly different 
(P<0.05) among treatments for that factor (or source of variability). 


Variability source df 


Onset hatching time 
S 
T 
Sx TE 
Error 
Total 
Corrected total 
Hatching rate 
S 
fp 
Sx T 
Error 
Total 
Corrected total 
Percentage of abnormality 
S 
T 
Ss x-r 
Error 
Total 
Corrected total 
Total length 
S 
ii) 
ox T 
Error 
Total 
Yolk sac volume 
S 
T 
SxT 
Error 
Total 
Survival rate 
S 
§ i 
ST 
Error 
Total 
Corrected total 


high mortality rates of late larvae and early juveniles due 
to impairments in prey capture and swimming behaviors. 

Yolk sac volume plays an important role during the 
endogenous feeding period that serves as the major source 
of energy for larvae, providing protein, free amino acids, 
lipids, and carbohydrates for maintenance, differentiation, 
growth, and other activities (Subhan et al., 2020). The 
results of our study indicate a significant downward trend 


Mean 
squares 


0.083 
LOR 
0.056 
0.057 


1.233 
333.251 
4.520 
95.907 


123.274 
249.836 
18.080 
82.840 


0.082 
0.015 
0.014 
0.004 


2.025 x 10° 
0.001 
1.582 x 107° 


208.333 
1075.000 
52.778 
190.000 


in yok sac volume of newly hatched larvae as turbidity 
increased, but no significant differences in volume were 
observed among salinity treatments or among treatments 
with different interactions between salinity and turbidity. 
Increased turbidity during the early developmental stages 
of embryos could decrease yolk sac volume, resulting 
in starvation of larvae before the mouth is fully formed 
for feeding and a high risk of mortality. Kamler (2008) 
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suggested that endogenous factors like yolk sac volume 
are one of the major causes of early mortality in fish. In 
addition, Subhan et al. (2020) reported that yolk sac vol- 
ume is a main factor in determining the quality of larvae 
in the postlarval stage that is critical to adaptation in the 
natural environment. 

Total length of newly hatched larvae is the only char- 
acteristic for which individual and combined effects of 
salinity and turbidity were found. Although stored energy 
is usually used for growth and organ maturation in nor- 
mal conditions, under stressful conditions, such energy 
is reallocated to cope with stress instead of being used 
for development (Pérez-Robles et al., 2016). In our study, 
total length of newly hatched larvae decreased when 
embryos were exposed to combinations of low salinity and 
high turbidity, indicating that embryos in non-optimal 
salinities may be unable to maintain their osmoregula- 
tory capacity because of decreasing availability of meta- 
bolic energy, thereby decreasing growth in development. 
Fish larvae with small sizes (<2.57 mm total length) at 
hatching have a lower chance of survival than larger 
larvae because they have fewer energy reserves, less 
feeding ability, and a slower, weaker response to preda- 
tors (Garrido et al., 2015). Furthermore, larvae that are 
smaller than other larvae at hatching cannot grow and 
survive because of a lower metabolic rate and higher cost 
of development. 


Conclusions 


Results from analysis of all characteristics of embryos 
and larvae measured in this study indicate that, although 
embryos of bamboo sole are quite tolerant and are able 
to withstand conditions of low salinity, they are sensitive 
to stress from high turbidity. As the frequency of heavy 
rainfall events continues to increase, reductions in the 
hatching rate of embryos and the size and survival rate of 
larvae of bamboo sole could lead to detrimental effects at 
the population level. 


Resumen 


Las lluvias intensas pueden disminuir la salinidad y 
aumentar la turbidez del agua en las zonas costeras, lo 
que afecta negativamente al desarrollo de los organis- 
mos, especialmente durante sus primeras etapas de vida. 
En este estudio, se evaluaron los efectos de la salinidad y 
la turbidez en los embriones del lenguado bambt (Hetero- 
mycteris japonicus) para comprender mejor su tolerancia 
al cambio climatico global. Se realizaron tres experimen- 
tos durante un periodo de 7 d. En el primer experimento, 
los embriones de lenguado bambt se expusieron durante 
3h a 1 de los 6 niveles de salinidad (14-84). Las longi- 
tudes totales de las larvas recién eclosionadas a niveles 
de salinidad bajos (14 y 18) fueron mas pequenas respecto 
al tamano de las larvas en los tratamientos con may- 
ores salinidades, aunque no se observaron diferencias 
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significativas en las tasas de eclosi6n y de supervivencia 
de las larvas entre los tratamientos. En el segundo exper- 
imento, los embriones se expusieron a turbideces de 0, 
100, 300, 500 y 700 unidades nefelométricas de turbidez 
durante 3 h. La turbidez disminuy6 significativamente 
la tasa de eclosién, la tasa de supervivencia y la longi- 
tud total, y aumento el tiempo de eclosion y el porcentaje 
de anormalidad. En el tercer experimento, los embriones 
se expusieron a diferentes combinaciones de salinidad 
y turbidez. El efecto de interaccién de la salinidad y 
la turbidez sobre la longitud total de las larvas recién 
eclosionadas fue significativo. Estos resultados indican 
que el desarrollo embrionario del lenguado bambt se vio 
mas afectado por los cambios de turbidez que por los de 
salinidad. 
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Guidelines for authors 


Contributions published in Fishery Bulletin describe 
original research in marine fishery science, fishery engi- 
neering and economics, and marine environmental and 
ecological sciences (including modeling). Preference will 
be given to manuscripts that examine processes and 
underlying patterns. Descriptive reports, surveys, and 
observational papers may occasionally be published but 
should appeal to an audience outside the locale in which 
the study was conducted. 

Although all contributions are subject to peer review, 
responsibility for the contents of papers rests upon the 
authors and not on the editor or publisher. Submission 
of an article implies that the article is original and is not 
being considered for publication elsewhere. 

All submissions are subject to a double-blind review 
to remove the identity of author and reviewer during the 
review process. 

Plagiarism and double publication are considered 
serious breaches of publication ethics. To verify the orig- 
inality of the research in papers and to identify possible 
previous publication, manuscripts may be screened with 
plagiarism-detection software. 

Manuscripts must be written in English; authors 
whose native language is not English are strongly advised 
to have their manuscripts checked by English-speaking 
colleagues before submission. 

Once a paper has been accepted for publication, online 
publication takes approximately 2 weeks. 

There is no cost for publication in Fishery Bulletin. 


Types of manuscripts accepted by the journal 


Articles generally range from 20 to 30 double-spaced 
typed pages (12-point font) and describe an original 
contribution to fisheries science, engineering, or econom- 
ics. Tables and figures are not included in this page count, 
but the number of figures should not exceed 1 figure for 
every 4 pages of text. Articles contain the following divi- 
sions: abstract, introduction, methods, results, discussion, 
and conclusions. 


Notes are generally less than 10 double-spaced typed 
pages (12-point font), including the Literature cited sec- 
tion. Like articles, notes describe an original contribu- 
tion to fisheries science. They follow a format similar to 
that for articles: abstract, introduction, methods, results, 
and discussion, but the results and discussion sections 
may be combined and a conclusions section should not 
be included. They should include no more than 2 figures 
or tables (2 of each would be too many). They are distin- 
guished from full articles in that they report a noteworthy 


new observation or discovery—such as the first report of 
a new species, a unique finding, condition, or event that 
expands our knowledge of fisheries science, engineering, 
or economics—and do not require a lengthy discussion. 
Manuscripts on range extensions will not be considered. 


Companion articles should be submitted together and 
are published together as a scientific contribution. Both 
articles should address a closely related topic and may be 
articles that result from a workshop or conference. 


Review articles of exceptional quality may be consid- 
ered occasionally for publication on a case-by-case basis. 
They address a timely topic that is relevant to aspects 
of fisheries science. They should include an abstract, but 
the format of the article, per se, will be up to the author. 
Please contact the scientific editor to discuss your ideas 
regarding a potential review article before embarking on 
such a project. 


Perspective essays generally are less than 5 double- 
spaced typed pages (12-point font) and provide clarifica- 
tion of misused terms or misunderstood topics in fisheries 
science. Essays are accepted at the discretion of the 
scientific editor. Opinion pieces on policy will not be con- 
sidered. They should include an abstract and a list of the 
literature cited, but no other sections should divide the 
article. Please contact the scientific editor to discuss your 
idea before making a submission. 


Preparation of manuscript 


Title page should include authors’ full names, affiliations, 
mailing addresses, and the senior author’s email address. 


Abstract should be limited to 200 words (one-half typed 
page), state the main scope of the research, and empha- 
size the authors’ conclusions and relevant findings. Do 
not review the methods of the study or list the contents of 
the paper. Because abstracts are circulated by abstract- 
ing agencies, it is important that they represent the 
research clearly and concisely. 


General text must be typed in 12-point Times New 
Roman font throughout. A brief introduction should 
convey the broad significance of the paper; the remain- 
der of the paper should be divided into the following 
sections: Materials and methods, Results, Discussion, 
Conclusions, and Acknowledgments. Headings within 
each section must be short, reflect a logical sequence, 
and follow the rules of subdivision (i.e., there can be no 
subdivision without at least 2 subheadings). The entire 
text should be intelligible to interdisciplinary readers; 
therefore, all acronyms, abbreviations, and technical 
terms should be written out in full and defined the first 
time they are mentioned. Abbreviations should be used 
sparingly because they are not carried over to index- 
ing databases and slow readability for those readers 
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outside a discipline. They should never be used for the 
main subject (e.g., species or method) of a paper. 

For general style, follow the U.S. Government Publish- 
ing Office Style Manual (2016, available at website) and 
Scientific Style and Format: the CSE Manual for Authors, 
Editors, and Publishers (2014, 8th ed.) published by the 
Council of Science Editors. For scientific nomenclature, 
use the current edition of the American Fisheries Soci- 
ety’s (AF'S) Common and Scientific Names of Fishes from 
the United States, Canada, and Mexico and its companion 
volumes (Crustaceans, Mollusks, Cnidaria and Ctenophora, 
and World Fishes Important to North Americans). For spe- 
cies not found in the previously mentioned AFS publica- 
tions and for more recent changes in nomenclature, use the 
Integrated Taxonomic Information System (ITIS, available 
at website), or, secondarily, the California Academy of Sci- 
ences Catalog of Fishes (available at website) for species 
names not included in ITIS. Common (vernacular) names of 
species should be lowercase. Citations must be given for the 
identification of specimens. For example, “Fish species were 
identified according to Collette and Klein-MacPhee (2002); 
sponges were identified according to Stone et al. (2011).” 

Dates should be written as follows: 11 November 2018. 
Measurements should be expressed in metric units, for 
example, “58 metric tons (t);” if other units of measurement 
are used, please make this fact explicit to the reader. Use 
numerals, not words, to express whole and decimal num- 
bers in the general text, tables, and figure captions (except 
at the beginning of a sentence). For example, “We consid- 
ered 3 hypotheses. We collected 7 samples in this location.” 
Use American spelling. Refrain from using the shorthand 
slash (/), an ambiguous symbol, in the general text. 

Cite all software, special equipment, and chemical solu- 
tions used in the study within parentheses in the general 
text, including the version number, company name, and the 
city and state (or nation) of the company headquarters, for 
example, “SAS, vers. 9.4 (SAS Institute Inc., Cary, NC).” 


Word usage and grammar that may be useful are the 
following: 


° Aging 
For our journal, the word aging is used to mean both 
age determination and the aging process (senescence). 
Authors should make clear which meaning is intended 
where ambiguity may arise. 


e Fish and fishes 
The plural of the word fish (a collective noun that 
implies individuals without regard to species) is fish. 
Example: The fish were collected by trawl net. 
Example: The numbers of fish collected that season 
were less than the numbers from previous years. 


The plural for fish species is fishes (a contrived plural 
used by taxonomists to mean several or more fish spe- 
cles) or one can use fish species (which is preferred in 
this journal for clarity across disciplines). 
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Example: The fishes of Puget Sound [biodiversity is 
implied] or 

Example: The fish species of Puget Sound [preferred 
plural for clarity across disciplines]. 


e Crab and crabs, squid and squids, etc. 
The plural of the word crab (i.e., many individuals 
without regard to species) is crab. 
Example: The crab were sorted by weight. 
Example: Many red king crab were dying [Many 
individuals of one species of crab.] 


The plural of crab species is crabs (a word used by tax- 
onomists) or crab species (the latter is preferred in this 
journal for clarity). 
Example: These crabs were selected for treatment. 
[Different crab species are implied.| 
Example: These crab species were selected for 
treatment. [Preferred word choice for clarity.] 
Example: Snow crabs are found throughout the 
North Pacific Ocean and Bering Sea. [There are 
2 species of snow crab; therefore the word crabs can 
be used here. | 
Example: Two species of snow crab are found through- 
out the North Pacific Ocean and Bering Sea. [Pre- 
ferred usage for clarity.] 
Example: Three crabs were selected for treatment. 
[3 species of crab are implied.] 
Example: Three crab species were selected for treat- 
ment. [Preferred word choice for clarity.] 


e We use fisherman and fishermen, not fisher and fishers, 
in this journal. One can always use crew member, vessel 
operator, and angler (the latter for recreational fishing). 


e The definite article with common names of species 
When the singular common name of a species rep- 
resents the entire class or group to which it belongs, 
use the definite article. 

Example: Only one species of the genus Salmo is 
found in the Atlantic Ocean—the Atlantic salmon 
(Salmo salar). 

Example: The sonic emissions of the bottlenose dol- 
phin are complex. 


For plural common names, this rule does not apply. 
Example: Chinook salmon are found throughout the 
Pacific Ocean. 
Example: Bottlenose dolphins are found in temper- 
ate and tropical waters. 


e Sex 
For the meaning of male and female, use the word sex, 
not gender. Do not write, “fish were sexed.” Write, “sex 
was determined.” 


e Participles 
As adjectives, participles must modify a specific noun 
or pronoun. 
Example: Using mark-recapture methods, these sci- 
entists determined the size of the population. [Correct. 
The participle wsing modifies the word scientists.] 
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Example: These scientists, based on the collected 
data, concluded that the mortality rate of these fish 
had increased. [Incorrect. The scientists were not 
based on the collected data.] 

Example: These scientists concluded, on the basis of 
collected data, that the mortality rate of these fish 
had increased. [Correct. The offending participle has 
been eliminated and an adverbial phrase modifies 
the verb concluded.| 


Equations and mathematical symbols should be set from a 
standard mathematical program (MathType or Equation 
Editor). Equations formatted in LaTex are not accept- 
able. For mathematical symbols in the general text (a, x’, 
T, +, etc.), use the symbols provided by the MS Word pro- 
gram and italicize all variables, except those variables 
represented by Greek letters and the superscript and 
subscript parts of variables and expressions. Do not use 
photo mode when creating these symbols in the general 
text, and do not cut and paste equations, letters, or sym- 
bols from a different software program. 

Number equations (if there are more than one) for 
future reference by scientists; place the number within 
parentheses at the end of the first line of the equation. 


Literature cited section comprises published works and 
those accepted for publication (in press) in peer-reviewed 
journals. Follow the name and year system for citation for- 
mat in this section (i.e., citations should be listed alpha- 
betically by the authors’ last names, and then by year if 
there is more than one citation by the same author. A list 
of abbreviations for citing journal titles can be found on 
our website. 

Authors are responsible for the accuracy and com- 
pleteness of all citations. Avoid the use of multiple 
citations when a single citation sufficiently supports a 
statement; cite the work that first reported the informa- 
tion that supports a statement, not all of the subsequent 
works. Literature citation format: Authors (last name, 
followed by initials for first name and, if given, mid- 
dle name of first author; then list names of additional 
authors with initials before last names). Year. Title of 
article. Abbreviated title of the journal in which it was 
published. Always include either the range of page num- 
bers (for a journal article) or a total number of pages (for 
a book or other type of publication). List a sequence of 
citations in the general text chronologically, for example, 
“(Smith, 1932; Green, 1998; Smith and Jones, 2015).” 


Acknowledgments should be no more than 6 lines of 
text. Only those who have contributed in an outstanding 
way should be acknowledged by name. For recognition of 
other persons or groups, use a general term, such as crew, 
observers, or research coordinators, and do not include 
names with these terms. 


Digital object identifier (doi) code ensures that a publica- 
tion has a permanent location online. A doi link (which 


may include a doi code) should be included at the end of 
citations of published literature. Authors are responsi- 
ble for submitting accurate doi links. Faulty links will be 
deleted at the page-proof stage. 


Footnotes are used for all documents that have not been 
formally peer reviewed and for observations and personal 
communications, but these types of references should be 
cited sparingly in manuscripts submitted to the journal. 

All reference documents, administrative reports, inter- 
nal reports, progress reports, project reports, contract 
reports, personal observations, personal communications, 
unpublished data, manuscripts in review, and council meet- 
ing notes are footnoted in 10-point font and placed at the 
bottom of the page on which they are first cited. Footnote 
format is the same as that for formal literature citations. A 
link to the online source (e.g., [Available from http://www... , 
accessed July 2017.]), or the mailing address of the agency 
or department holding the document, should be provided 
so that readers may obtain a copy of the document. 


Tables are often overused in scientific papers; it is seldom 
necessary to present all the data associated with a study. 
Tables should not be excessive in size and must be cited 
in numerical order in the text. Headings should be short 
but ample enough to allow the table to be intelligible on 
its own. 

All abbreviations and unusual symbols must be 
explained in the table legend. Other incidental com- 
ments may be footnoted with numeral footnote markers. 
Use asterisks only to indicate significance in statistical 
data. Do not put a table legend on a page separate from 
the table; place the legend above the table. Do not submit 
tables in photo mode. 


e Note probability with a capital, italic P. 


e Provide a zero before all decimal points for values less 
than one (e.g., 0.07). 


e Round all values to 2 decimal points. 


e Use a comma in numbers of 5 digits or more (e.g., 
13,000 but 3000). 


Figures must be cited in numerical order in the text. 
Graphics should aid in the comprehension of the text, but 
they should be limited to presenting patterns rather than 
raw data. The number of figures should not exceed 1 figure 
for every 4 pages of text. 

Figure legends should explain all symbols and abbrevi- 
ations seen in the figure and should be double spaced on a 
separate page at the end of the manuscript. 

Line art and halftone figures should be saved at res- 
olutions >600 dots per inch (dpi) and >300 dpi, respec- 
tively. Color is allowed in figures to show morphological 
differences among species (for species identification), to 
show stain reactions, to show gradations (such as those of 
temperature and salinity within maps), and to distinguish 
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between numerous lines and symbols in graphs. Figures 
approved for color should be saved in CMYK format. 

All figures must be submitted as PDF, TIFF, or EPS 
files. 3 


e Capitalize only the first letter of the first word and 
proper nouns in all labels within figures. 


e Do not use overly large font sizes for labels in maps and 
for axis labels in graphs. 


e Use the same point size for all labels, except for panel 
labels (e.g., A), which should be slightly larger than 
other labels (e.g., 11 point versus 8 point). 


e Use a sans serif font for all labels. Panel labels (e.g., A), 
however, should be in Times New Roman font. 


¢ Do not use bold fonts or bold lines in figures. Do not 
use italic fonts. Exceptions include use of italic fonts 
for labels of bodies of water in maps and a bold font for 
panel labels (e.g., A). 


e Do not place outline rules around graphs. 


e Do not include vertical and horizontal lines in the back- 
ground of graphs. Ticks for values on the x-axis and 
y-axis will suffice. 


e Place a north arrow and label degrees latitude and lon- 
gitude (e.g., 170°E) in all maps. If scale of map requires 
more than degrees, use degrees minutes, not decimal 
degrees. 


e Place panel labels (e.g., A, B, C) within the upper-left 
area of each graph or photo in a multi-panel figure, 
from left to right, then top to bottom. If the letter is 
not visible against a dark background, put a white box 
behind it. Do not use white labels. 


e Avoid placing labels vertically or diagonally. Y-axis 
labels can be vertical. Words in horizontal labels can be 
stacked vertically to fit. 


e Use symbols, shading, or patterns (not clip art) in maps 
and graphs. 


e For scale bars in maps, use kilometers. Use the label 
km or kilometers (lowercased). 


Supplementary materials that are considered essential, 
but are too large or impractical for inclusion in a paper 
(e.g., metadata, figures, tables, videos, and websites), may 
be provided at the end of an article. These materials are 
subject to the editorial standards of the journal. A URL to 
the supplementary material and a brief explanation for 
including such material should be sent at the time of ini- 
tial submission of the paper to the journal. 


e Metadata, figures, and tables should be submitted in 
standard digital format (MS Word or PDF file) and 
should be cited in the general text, for example, as 
“.. was determined (Suppl. Table 3, Suppl. Fig. 1).” 


e Websites should be cited with a URL in the general 
text. 


e Videos must not be larger than 30 MB to allow a swift 
technical response for viewing the video. Authors 
should consider whether a short video uniquely cap- 
tures what text alone cannot capture for the under- 
standing of a process or behavior under examination 
in the article. Supply an online link to the location of 
the video. 


Copyright law does not apply to Fishery Bulletin, which 
falls within the public domain. However, if an author 
reproduces any part of an article from Fishery Bulletin, 
reference to source is considered correct form (e.g., Source: 
Fish. Bull. 117:105). 


Failure to follow these guidelines 
and failure to correspond with editors 
in a timely manner will delay 
publication of a manuscript. 


Submission of manuscript 


Submit a manuscript online at the ScholarOne Manu- 
scripts website for Fishery Bulletin. Commerce Depart- 
ment authors must provide proof of internal clearance 
of their manuscripts with either a completed and signed 
NOAA Form 25-700 or a copy of the clearance email from 
the Research Publication Tracking System. For further 
details on electronic submission, please contact the asso- 
ciate editor, Cara Mayo, at 


cara.mayo@noaa. gov. 


When requested, the text and tables should be submitted 
in MS Word format. Each figure should be sent as a sep- 
arate PDF, TIFF, or EPS file. Send a copy of a figure in 
the original software if conversion to any of these formats 
yields a degraded version of the figure. 


Questions? If you have questions regarding these guide- 
lines, please contact the managing editor, Kathryn 
Dennis, at 


kathryn.dennis@noaa. gov. 


Questions regarding manuscripts under review should be 
addressed to Associate Editor Cara Mayo. 
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